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Abstract 


The  purpose  of  this  research  was  to  describe  the  unperturbed  relative  motion  of  Earth  satellites  in 
elliptical  orbits  using  a  simple  dynamics  model  whose  parameters  allow  significant  geometrical  insight  and 
operational  efficacy.  The  goal  was  to  retain  the  advantages  of  the  Relative  Orbit  Elements  (ROE) 
realization  of  the  Hill-Clohessy-Wiltshire  (HCW)  equations,  a  linearized  dynamics  model  for  circular 
reference  orbits.  Specifically,  this  thesis  analyzed  the  geometry  of  satellite  rendezvous  and  proximity 
operations  using  the  ROE  parameters  to  characterize  the  model’s  utility.  Next,  through  a  comprehensive 
literature  review,  this  thesis  sought  possible  approaches  for  developing  a  similarly  useful  parameterization 
for  chief  orbits  with  nonzero  eccentricity.  The  approach  selected  was  a  novel  linear  time-varying  system 
which  requires  both  chief  and  deputy  satellites  to  remain  close  to  a  virtual  chief  on  a  known  circular  orbit. 
The  research  derived  and  solved  the  equations  of  motion,  expressing  the  solution  in  terms  of  simple 
geometric  parameters.  Numerical  simulations  compared  the  new  model  against  both  HCW  and  Keplerian 
two-body  motion,  revealing  less  accurate  performance  than  HCW  for  some  cases.  Error  analysis  explained 
this  behavior  and  found  restricted  regions 
identified  new  approaches  for  researching 


where  the  new  model  performed  accurately.  Einally,  this  study 
relative  satellite  motion  on  elliptical  orbits. 
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RELATIVE  ORBIT  ELEMENTS  FOR  SATELLITES  IN  ELLIPTICAL  ORBITS 


I.  Introduction 

Background;  Relative  Satellite  Motion  Models 

Dynamics  models  are  mathematical  representations  of  how  an  object  or  system  of 
objects  behaves  in  accordance  with  the  physics  of  forces  and  motion.  A  relative  satellite 
motion  model  is  a  type  of  dynamics  model  describing  the  motion  of  one  orbiting 
spacecraft  as  seen  from  another.  Modeling  the  relative  motion  of  Earth  satellites  is  a 
tremendously  important  subject  with  a  great  diversity  of  applications,  as  we  shall  see.  It 
is  also  a  well-studied  subject;  many  important  contributions  to  the  study  were  made 
before  and  during  the  dawn  of  the  Space  Age. 

In  fact,  the  most  celebrated  accomplishments  of  the  Space  Age  would  not  have 
been  possible  without  a  good  understanding  of  relative  satellite  motion.  The  Apollo 
astronauts  could  not  have  traveled  safely  to  the  Moon  and  back  without  the  capability  to 
rendezvous  and  dock  in  their  spacecraft;  this  is  a  key  application  of  relative  satellite 
motion.  Likewise,  the  history  of  international  cooperation  in  space,  dating  back  to  the 
Apollo-Soyuz  missions,  has  depended  to  some  extent  on  merging  two  very  different 
technologies  for  achieving  spacecraft  rendezvous  and  docking:  the  astronaut-dependent 
approach  favored  by  the  United  States  and  the  more  automated  approach  favored  by 
Russia  (and,  formerly,  the  Soviet  Union). 
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Today’s  space  technology  involves  many  other  practical  uses  of  relative  satellite 
motion.  These  include  formation  flying,  in  which  multiple  spacecraft  orbit  together  with 
their  relative  motion  precisely  controlled  in  order  to  obtain  and  correlate  certain  data  from 
the  environment.  They  also  include  proximity  operations,  in  which  one  satellite 
maneuvers  near  another,  possibly  for  purposes  other  than  docking,  such  as  photographing 
or  inspecting  for  damage  from  various  perspectives.  Proximity  operations  are  often 
grouped  together  with  techniques  for  terminal  rendezvous  in  a  category  of  operations 
called  RPO  (rendezvous  and  proximity  operations).  Although  RPO  of  various  types  have 
been  performed  for  decades,  it  has  only  been  recently,  and  with  considerable  difficulty, 
that  RPO  have  been  attempted  autonomously,  meaning  that  the  maneuvers  are  planned  by 
computers  on-board  one  or  more  of  the  spacecraft,  rather  than  by  humans  or  higher- 
powered  computers  on  the  ground.  As  we  shall  see,  autonomy  adds  special  requirements 
to  the  type  of  relative  satellite  motion  model  which  must  be  employed. 

At  this  point,  it  may  be  helpful  to  emphasize  the  subtle  but  important  distinction 
between  the  terms  automated  and  autonomous,  as  used  in  this  study.  Autonomy 
(sometimes  called  full  autonomy  or  autonomous  planning)  implies  that  guidance, 
navigation,  and  control  are  accomplished  with  neither  astronauts  nor  ground  stations  in 
the  loop,  and  that  maneuvers  are  planned  by  the  on-board  computer  [38,  44,  104]. 
Automated  (or  automatic)  RPO  has  been  used  since  1967,  when  early  unmanned  Soyuz 
vehicles  docked  during  the  Kosmos  186/188  mission  [104].  Automated  rendezvous  as 
still  practiced  by  Progress  cargo  vehicles  involves  a  pre-determined  flight  profile,  with 
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some  maneuvers  initiated  by  ground  controllers,  and  other  maneuvers  calculated  from  the 
relative  navigation  sensor  data  [31]. 

Selecting  the  best  dynamics  model  to  use  for  a  particular  application  often 
depends  on  a  trade-off  between  simplicity  and  accuracy.  Every  dynamics  model  is  an 
approximation;  that  is,  mathematically  representing  every  force  and  the  motion  of  every 
object  in  the  universe  is  neither  possible  nor  helpful.  Rather,  an  engineer  need  only 
model  those  forces  which  make  a  noticeable  difference  for  the  given  application. 
Likewise,  an  engineer  can  represent  those  forces,  as  well  as  the  resulting  motion,  with  a 
mathematical  approximation.  This  approximation  may  be  simpler  and  less  accurate  than 
the  most  accurate  description  offered  by  physics,  so  long  as  the  resulting  error  is  not  too 
great  for  the  given  application. 

As  an  example,  consider  a  spacecraft  attempting  autonomous  RPO.  Autonomy  in 
space  systems  depends  on  the  computer  software,  including  navigation  filters,  guidance 
algorithms,  and  control  systems,  operating  without  direction  from  control  stations  on 
Earth.  On  the  one  hand,  the  dynamics  model  at  the  core  of  each  of  these  software 
functions  must  be  accurate  enough  for  any  given  operation  to  succeed  within  the 
mission's  error  tolerance.  On  the  other  hand,  if  the  model  is  not  simple  enough  to 
implement  with  the  available  computing  power,  the  spacecraft  cannot  perform  the 
function.  In  particular,  spacecraft  trajectory  design  is  usually  the  subject  of  extensive  pre¬ 
mission  planning  using  high-speed  computers  and  high-fidelity  physics  models.  This 
practice  is  good  and  useful,  but  for  truly  autonomous  RPO,  the  trajectory  planning 
function  must  be  replicated  on-board. 
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The  field  of  relative  satellite  motion  includes  models  that  consider  a  variety  of 
forces.  In  the  present  study,  though,  we  will  focus  on  models  which  include  only  a 
simplified  representation  of  the  Earth’s  gravity.  Each  satellite  will  be  considered  a  point 
mass,  and  the  Earth  will  be  considered  a  uniform  sphere,  so  that  its  gravity  field  decreases 
with  the  distance  squared  from  the  Earth’s  center.  We  call  this  approximation  the  two- 
body  assumption,  because  it  is  the  basis  for  the  classical  two-body  problem  of 
astrodynamics.  The  trajectory  of  a  spacecraft  moving  in  such  an  inverse-square  gravity 
field  is  a  simple  Keplerian  orbit,  which  can  be  described  in  terms  of  the  classical 
parameters  known  as  orbit  elements.  Commonly  known  properties  of  two-body  orbits 
can  be  found  in  any  astrodynamics  text;  for  convenience,  several  are  listed  in  Appendix 
D. 

Of  course,  in  reality,  satellites  experience  many  forces  other  than  inverse-square 
gravity.  In  general,  the  accuracy  of  the  two-body  assumption  tends  to  decrease  over  time 
as  these  forces  perturb  the  orbit.  However,  relative  satellite  motion  models  which  do  not 
consider  such  perturbing  forces  are  accurate  enough  for  some  applications  over  limited 
time  scales. 

Other  approximations  can  be  used  to  simplify  relative  satellite  motion  models. 
Many  models  assume  that  the  dynamics  are  governed  by  a  set  of  linear  differential 
equations;  this  linearization  usually  requires  assuming  small  differences  in  the  variables. 
Eor  example,  if  the  distance  between  two  satellites  is  small  enough  compared  to  their 
distance  from  the  Earth’s  center,  or  if  the  differences  in  their  classical  orbit  elements  are 
small  enough,  then  any  terms  in  the  equations  of  motion  which  contain  one  small  variable 
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multiplied  by  another  may  be  neglected.  Such  linearized  models  are  the  main  emphasis 
of  the  present  study.  If  a  given  application  cannot  accept  the  error  that  results  from 
neglecting  all  nonlinear  terms,  it  may  require  a  model  which  retains  the  quadratic  terms, 
but  neglects  terms  which  are  third-order  or  higher  in  the  variables;  these  quadratic  models 
are  more  accurate  but  also  more  complex  than  linearized  ones. 

Another  approximation  commonly  used  for  relative  satellite  motion  models 
involves  the  eccentricity  of  the  satellite  orbits.  Several  properties  of  two-body  orbits  can 
be  approximated  by  power  series  expansions  in  terms  of  orbital  eccentricity.  If  a  given 
application  can  be  restricted  to  satellites  in  low -eccentricity  orbits,  then  truncating  these 
power  series  expansions  after  the  first-order  or  second-order  terms  will  introduce  only  a 
small  amount  of  error.  If  the  application  is  restricted  to  orbits  with  negligible 
eccentricity,  then  only  the  zeroth-order  terms  need  to  be  retained. 

The  simplest  and  most  commonly-applied  relative  satellite  motion  models  address 
the  circular  chief  problem;  that  is,  they  describe  the  motion  of  one  satellite,  a  deputy,  with 
respect  to  a  chief  satellite  whose  orbit  has  zero  eccentricity.  Modeling  the  circular  chief 
problem  is  discussed  in  more  detail  in  Chapters  n  and  HI.  For  now,  it  suffices  to  say  that 
the  most  common  model  is  the  Hill-Clohessy-Wiltshire  (HCW)  model,  a  two-body 
linearized  model  [40,  16].  One  particularly  useful  realization  of  the  HCW  model  is  the 
Relative  Orbit  Elements  (ROEs)  developed  at  the  Air  Eorce  Research  Eaboratory  and 
AEIT  [60] .  As  we  shall  see,  ROEs  provide  unique  insight  into  the  geometry  of  the 
relative  trajectory  and  are  simple  enough  to  serve  as  the  foundation  for  maneuver 
schemes  and  guidance  algorithms  suitable  for  autonomous  RPO.  Eovell  and  Tragesser 
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described  this  advantage  of  ROEs  as  "operational  efficacy",  because  they  "lend 
themselves  well  to  the  development  of  on-board  guidance  and/or  mission  planning 
software"  [60]. 

Another  group  of  relative  satellite  motion  models  address  the  elliptical  chief 
problem,  in  which  the  chiefs  orbit  is  permitted  to  be  noncircular;  that  is,  the  chief 
eccentricity  may  take  on  values  between  zero  and  one.  Modeling  the  elliptical  chief 
problem  is  discussed  in  more  detail  in  Chapters  n  and  IV.  In  this  problem,  even  after 
making  the  two-body  assumption,  a  satellite’s  angular  velocity,  angular  acceleration,  and 
distance  from  the  Earth’s  center  are  no  longer  constants,  as  they  are  in  the  circular  chief 
problem.  Eor  this  reason,  elliptical  chief  models  tend  to  be  more  complex. 

Problem  Statement 

Although  many  relative  satellite  motion  models  exist  for  the  elliptical  chief 
problem,  none  have  been  stated  in  a  way  that  lends  itself  to  the  same  high  degree  of 
operational  efficacy  as  have  ROEs  for  the  circular  chief  problem.  In  other  words,  there  is 
a  need  for  realizations  of  relative  satellite  motion  which  do  not  restrict  the  eccentricity  of 
the  chief  satellite’s  orbit  to  zero,  but  which  are  still  practical  enough  to  facilitate 
autonomous  guidance  and  maneuver  schemes. 

This  need  is  the  basis  for  the  current  research.  This  thesis  attempts  to  answer  the 
following  two  questions:  What  characteristics  of  the  ROE  realization  provide  operational 
efficacy?  How  can  we  retain  those  advantageous  characteristics  while  relaxing  the 
eccentricity  restriction? 
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The  thesis  is  organized  in  five  chapters.  The  remainder  of  Chapter  I  introduces 
the  reference  frames  and  some  of  the  mathematical  notation  used  throughout  the  later 
chapters.  Chapter  n  is  a  comprehensive  (although  not  exhaustive)  survey  of  the  vast 
literature  on  relative  satellite  motion  models.  Chapter  HI  includes  a  brief  overview  of  the 
HCW  model  and  a  detailed  derivation  and  geometric  analysis  of  the  ROE  parameters  for 
the  circular  chief  problem.  Chapter  IV  lays  out  five  possible  approaches  for  handling 
eccentricity  in  the  chiefs  orbit;  it  then  follows  through  with  the  approach  selected,  called 
the  Virtual  Chief  method,  by  deriving  and  solving  the  equations  of  motion  for  a  new 
relative  satellite  motion  model.  The  chapter  continues  by  deriving  a  parameterized 
version  of  the  Virtual  Chief  model’s  solution,  analyzing  the  new  parameters  for  their 
geometric  significance.  Chapter  V  examines  numerical  simulation  results  and  evaluates 
the  accuracy  of  the  Virtual  Chief  model;  finding  its  accuracy  significantly  less  than  would 
be  desired,  the  chapter  concludes  with  a  theoretical  analysis  of  the  approximation  errors 
and  recommendations  for  future  development. 

Notation  and  Reference  Frames 

The  reference  frame  most  commonly  used  to  describe  relative  satellite  motion  is  a 
local-vertical  local-horizontal  (LVLH)  frame,  sometimes  called  an  orbital  or  Hill  frame. 
The  center  of  an  LVLH  frame  is  fixed  to  a  satellite,  and  its  coordinate  axes  rotate  at  the 
same  rate  as  the  satellite  moves  through  its  orbit.  When  considering  a  relative  satellite 
motion  model,  we  can  say  that  the  chief  satellite  exists  at  point  C,  and  its  LVLH  frame  is 
frame  c,  as  in  Ligure  1.  The  frame’s  orthogonal  unit  basis  vectors  are  defined  in  this 
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thesis  as  follows:  Cj  is  parallel  to  the  position  veetor  from  the  Earth's  center  to  the  chief, 
in  the  outward  sense.  Cj  is  orthogonal  to  ,  within  the  plane  of  the  chiefs  orbit,  and  in 
the  same  sense  as  the  chiefs  inertial  velocity  vector,  completes  the  right-handed  set. 


D^ 


Of  course,  dynamics  cannot  be  performed  without  an  inertial  reference.  In 
studying  relative  satellite  motion,  the  most  commonly  used  inertial  frames  are  Earth- 
centered  inertial  (ECI)  frames.  Some  of  these  are  perifocal  frames,  oriented  to  a 
satellite’s  orbit  plane.  If  i  is  an  ECI  frame,  then  we  can  say  that  the  angular  velocity  of 

frame  c  with  respect  to  i  is  a, . 

Yi 

Also,  if  the  deputy  satellite  is  located  at  point  D,  we  can  say  that  its  position 

relative  to  the  chief  is  .  A  relative  satellite  motion  model  must  be  able  to  describe  not 

/c 
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only  position,  but  also  the  deputy’s  relative  veloeity,  or  the  time  derivative  of  the  relative 


position  vector  with  respect  to  the  c-frame,  denoted  —  .  This  description  is  usually 

dt  /c 

(although  not  always)  accomplished  via  equations  of  motion  which  define  the  relative 

c  j2 

acceleration,  the  c-frame  second  time  derivative,  — ^  . 

dt  /c 

If  we  assign  variables  a: ,  y  ,  and  z  as  the  c-frame  relative  position  components, 
then  we  can  write  the  relative  position  in  either  of  the  following  ways: 

=  xcj  -1-  yc2  -f  zcj 
X 

=  y 

z 


The  solution  of  a  relative  satellite  motion  model,  then,  must  give  values  for  all 
three  components  of  relative  position  and  all  three  components  of  relative  velocity  for  any 
time  t .  In  other  words,  using  a  dot  to  represent  a  scalar  time  derivative,  we  must  know 
x(t),y(t)  ,  z[t)  x(t),  y(t),and  z(t).  Note  that 


X 


y 


z 


Ic 


The  coordinates  just  described  constitute  a  rectangular  Cartesian  coordinate 
system.  It  should  be  noted  that  cylindrical  coordinates  may  also  be  used  for  an  LVLH 


9 


frame;  in  this  case,  y  is  replaced  by  r^A.0 ,  where  is  the  distance  from  the  Earth’s 

center  to  the  chief,  and  A0  is  the  polar  angle  between  the  deputy  and  the  chief,  as  shown 
in  Figure  1.  Interpreting  an  LVLH  solution  as  cylindrical  can  increase  the  accuracy,  since 
the  chief  satellite  is  not  really  moving  in  a  straight  line  down  the  y-axis.  This  is 
especially  true  in  the  circular  chief  problem,  when  is  constant  and  9  corresponds  to 
the  true  longitude  through  which  the  chief  is  moving  at  a  constant  rate.  However,  the 
error  due  to  using  rectangular  coordinates  is  often  well  within  a  given  application’s  error 
tolerance,  particularly  if  the  application  is  using  a  linearized  model,  and  thus  already 
assuming  small  relative  separation  distances. 
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II.  Literature  Review 


There  is  a  vast  amount  of  technical  literature  on  relative  satellite  motion.  This 
survey  cannot  be  exhaustive,  but  it  is  intended  to  cover  the  most  significant  developments 
in  a  variety  of  areas.  The  first  two  sections  describe  models  of  unperturbed  motion  in  the 
circular  and  elliptical  chief  problems.  The  third  section  explains  some  of  the  most 
significant  practical  applications  that  have  been  designed  to  use  such  models.  The  final 
two  sections  provide  a  brief  survey  of  higher-fidelity  models  and  review  a  few  existing 
methods  for  making  comparisons  between  models. 

Circular  Chief  Problem 

Hill-  Clohessy- Wiltshire  Model. 

The  first  and  best  known  relative  satellite  motion  model  is  the  HCW  model 
mentioned  in  Chapter  I.  It  was  developed  by  Clohessy  and  Wiltshire  in  1960  [16].  The 
equations  of  motion  which  they  derived  and  solved  are  similar  to  those  used  by  the 
famous  American  mathematician  George  Hill  in  1878  to  describe  the  Moon's  motion 
relative  to  the  Earth  [40].  The  HCW  model  assumes  that  each  satellite  is  following 
Keplerian  two-body  motion  and  that  the  chiefs  orbit  is  circular;  the  equations  of  motion 
are  linearized  by  assuming  small  relative  distances.  The  HCW  model  has  served  as  the 
basis  for  most  real-life  RPO  missions.  As  we  shall  see,  a  key  advantage  of  HCW  is  that 
its  solution  has  easily  discernible  geometric  features. 

In  2001,  Sabol,  Bums,  and  McLaughlin  published  a  study  closely  examining  the 
geometry  of  HCW  dynamics  [75].  They  noted  that  even  very  small  differences  in  semi- 
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major  axis  between  the  ehief  and  deputy  would  result  in  different  periods  (see  Eq.  (D.l)) 
and  therefore  would  eause  the  two  satellites  to  drift  apart  in  the  y-direetion  over  time. 
They  further  asserted  that,  under  HCW  dynamies,  every  relative  satellite  trajeetory  whieh 
does  not  drift  in  this  way  is  shaped  like  an  ellipse.  They  also  proposed  a  number  of 
potentially  useful  formation  flying  geometries,  ineluding  a  projeeted  eireular  orbit  (PCO), 
in  whieh  several  satellites  orbit  together  in  sueh  a  way  that  they  appear  to  remain  on  a 
eirele  in  the  y-z  plane. 

In  2004,  Lovell,  Tragesser,  and  Tollefson  [61]  and  Lovell  and  Tragesser  [60] 
identified  the  ROEs,  six  parameters  that  eould  be  used  to  define  HCW  relative  satellite 
motion  in  a  geometrieally  intuitive  way.  They  showed  the  transformations  between  ROEs 
and  other  possible  realizations:  Cartesian  initial  eonditions  (that  is,  values  of  the  position 
and  veloeity  eomponents  at  epoeh  time)  and  orbit-element  differenees  (that  is,  differenees 
in  the  values  of  the  ehief  s  and  deputy's  elassieal  orbit  elements).  They  also  showed  that 
the  HCW  solution  written  in  terms  of  ROEs  ean  be  used  to  ealeulate  impulsive 
maneuvers  as  part  of  an  autonomous  guidanee  algorithm. 

Also  in  2004,  Press  used  ROEs  to  deseribe  the  physieal  geometry  of  an  HCW 
relative  orbit  [71].  He  showed  that  an  instantaneous  relative  trajeetory  is  eonfined  to  a 
plane  and  is  defined  by  the  interseetion  of  that  plane  with  a  2x1  elliptieal  eylinder 
orthogonal  to  the  x-y  plane.  The  relative  orbit  plane  and  the  elliptieal  eylinder  may  be 
drifting  at  the  same  rate  in  the  y-direetion. 
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Quadratic  Models. 

In  1963,  London  published  a  new  relative  satellite  motion  model  accurate  over 
longer  distances  than  the  HCW  model  [57].  He  still  assumed  a  circular  chief  orbit,  but 
retained  second-order  terms  in  the  relative  distance.  London's  equations  of  motion  are 
more  complex  than  the  HCW  equations,  but  they  can  be  solved  by  the  method  of 
successive  approximations.  The  resulting  solution  includes  secular  terms  (containing  t 
and  t^)m  the  y-direction  and  mixed-secular  terms  (of  the  form  a^t  sin  ( ajt  -i-  <23 ) )  in  all 
three  directions. 

In  1999,  Sedwick,  Miller,  and  Kong,  in  the  process  of  studying  HCW  dynamics 
perturbed  by  the  J2  term  of  the  Earth's  gravitational  field,  developed  an  unperturbed 
quadratic  equation  of  motion  for  the  radial  direction  [81].  They  illustrated  that,  for  some 
relative  satellite  motion  applications,  the  linearization  error  is  negligible  compared  to  the 
differential  J2  effect. 

More  recently,  in  2003,  Karlgaard  and  Lutze  used  spherical  instead  of  rectangular 
or  cylindrical  coordinates  to  derive  equations  of  relative  motion  [48].  They  retained  up  to 
quadratic  terms  in  the  relative  position  and  velocity  components.  They  solved  the 
equations  using  the  method  of  multiple  scales;  part  of  their  rationale  for  this  method  was 
that  it  eliminates  the  mixed-secular  terms  in  the  x-  and  z-directions  found  in  London's 
solution. 

About  the  same  time,  Gurfil  and  Kasdin  used  orbit-element  differences  to  create 
arbitrarily  high-order  nonlinear  corrections  to  the  HCW  model  [36,  37].  Then,  by 
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truncating  the  infinite  series  relating  the  deputy's  true  anomaly  to  time,  they  created  a 
time-explicit  solution  that  is  second-order  in  the  orbit-element  differences. 

Elliptical  Chief  Problem 

Lawden  's  Equations:  the  First  Generation. 

The  first  generation  of  linearized  models  for  the  elliptical  chief  problem  was 
developed  via  three  independent  paths  during  the  period  from  the  1950s  through  the  mid- 
1980s.  Only  one  of  the  three  paths  was  pursuing  a  relative  satellite  motion  model.  The 
other  two  paths  used  equivalent  equations  of  motion  while  attempting  to  model  other 
spaceflight  phenomena,  such  as  a  spacecraft's  deviation  from  a  reference  trajectory.  This 
system  of  equivalent  equations  has  been  expressed  in  a  variety  of  coordinate  types 
(rectangular  and  cylindrical)  and  independent  variables  (time,  chief  true  anomaly,  and 
chief  eccentric  anomaly).  Derivations  of  several  forms  of  these  equations  can  be  found  in 
Schaub  and  Junkins  [79]. 

Path  One. 

The  first  to  derive  and  solve  the  equations  was  Lawden  in  1954  (in-plane 
components  only)  [55]  and  1963  (all  three  components)  [56].  His  purpose  was  to 
describe  the  primer  vector  for  optimal  spacecraft  trajectory  design.  Lawden's  equations 
are  perhaps  the  best-known  elliptical-chief  model.  Written  for  LVLH-frame  coordinates, 
they  use  the  chiefs  true  anomaly  as  the  independent  variable,  so  that  the  first  equation 
defines  x" ,  a  scalar  second  derivative  with  respect  to  chief  true  anomaly,  rather  than  x  , 
and  so  on.  Lawden's  1963  book  lays  out  a  complete  theory  of  optimal  trajectory  design; 
the  sections  of  present  interest  concern  spacecraft  in  an  inverse-square  gravity  field. 
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Lawden  presented  a  closed- form  solution  for  position  and  velocity  in  terms  of  arbitrary 
integration  constants.  He  also  presented  the  equations  and  solution  in  an  alternate  form 
where  each  coordinate  is  scaled  by  the  chiefs  orbit  radius;  note  that  this  scale  factor  is  not 
a  constant  (see  Eq.  (D.2)).  When  the  solution  is  written  for  unsealed  coordinates,  its 
terms  contain  a  cumbersome  definite  integral  which  is  singular  for  circular  chief  orbits. 
However,  the  1954  article  gave  a  closed-form  expression  for  this  integral  that  is  valid  for 
chief  eccentricities  between  zero  and  one. 

A  key  characteristic  of  all  exact  solutions  to  Lawden's  equations  is  that  their 
coefficients  are  time-varying,  but  not  explicitly.  Instead,  the  coefficients  vary  with  the 
chiefs  true  anomaly.  Thus,  because  Kepler's  equation  (Eq.  (D.4))  is  transcendental,  the 
solutions  cannot  be  completely  explicit  in  time. 

A  little-noticed  contribution  from  Eawden's  book  was  his  study  of  the  primer 
vector  solution  in  coordinates  parallel  to  the  reference  perifocal  frame.  This  solution 
(which  is  mathematically  equivalent  to  studying  the  deputy's  relative  motion  in  a  non¬ 
rotating  frame  fixed  to  the  chief  and  parallel  to  the  chiefs  perifocal  frame)  proves  to  be 
much  simpler,  showing  a  small  variation  from  circular  motion. 

Path  Two. 

The  second  path  began  with  De  Vries,  also  in  1963  [26].  De  Vries  was  the  first  to 
address  the  elliptical  chief  problem  directly.  Eike  Clohessy  and  Wiltshire,  he  linearized 
the  two-body  problem  by  assuming  small  relative  distances  to  create  time-domain 
differential  equations  in  the  chiefs  EVEH  coordinates.  However,  because  he  did  not 
assume  a  circular  chief  orbit,  the  coefficients  of  these  equations  of  motion  include  time- 
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varying  quantities  (chief  orbit  radius  and  chief  true  anomaly).  In  order  to  solve  this 
system  of  equations,  De  Vries  attempted  to  transform  them  into  the  true-anomaly  domain. 
However,  the  out-of-plane  component  of  this  transformed  system  is  incorrect  [94] . 

Two  years  later,  Tschauner  and  Hempel  tackled  the  relative  satellite  motion 
problem  in  a  paper  whose  title  translates  to  "Rendezvous  with  a  Target  in  an  Elliptical 
Orbit"  [94].  They  corrected  De  Vries's  mistake  to  create  a  three-dimensional  system  of 
true-anomaly-domain  equations  of  relative  motion.  After  they  scaled  these  equations  to 
the  chiefs  orbit  radius,  they  could  be  solved  simply  in  closed  form.  Tschauner  and 
Hempel's  equations  of  motion,  in  both  scaled  and  unsealed  form,  are  mathematically  the 
same  as  those  used  by  Lawden.  In  fact,  some  authors  refer  to  Lawden's  equations  as  the 
Tschauner-Hempel  equations.  Later,  Tschauner  studied  the  equivalent  system  of 
equations  transformed  to  the  chief  eccentric  anomaly  domain  [93]. 

In  1966,  Shulman  and  Scott  presented  an  alternate  derivation  for  the  linearized 
equations  of  motion  [88].  They  proceeded  to  find  a  closed-form  solution.  They  also 
mapped  the  arbitrary  integration  constants  of  that  solution  to  Cartesian  initial  conditions 
for  the  special  case  when  epoch  time  occurs  at  chief  perigee.  Importantly,  they  published 
graphical  results  of  numerical  simulations  comparing  their  model  (a  linearized  model 
complete  in  the  chiefs  eccentricity,  like  Lawden's  equations)  to  the  full  nonlinear 
Keplerian  equations  integrated  numerically  and  to  a  variety  of  low -eccentricity 
approximations.  Some  of  these  numerical  results  show  the  effect  of  the  mixed-secular 
terms  in  the  x-direction. 
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In  1983,  Wolfsberger,  WeiB,  and  Rangnitt  presented  a  general  closed-form 
solution  to  Lawden's  equations  (which  they  knew  only  from  Tschauner  and  Hempel) 
[105].  These  solutions  take  the  form  of  a  linear  system  in  the  Cartesian  initial  conditions, 
and  can  be  initialized  with  an  epoch  time  anywhere  along  the  chief  orbit.  They  continue 
to  use  rotating  LVLH  coordinates  primarily,  although  they  did  calculate  and  plot  an 
example  relative  trajectory  in  inertial  chief-fixed  coordinates,  allowing  a  visual 
comparison. 

Path  Three. 

The  third  path  began  with  Stem's  research  in  1963  [90].  Stem  independently 
developed  a  model  mathematically  equivalent  to  De  Vries's  time-domain  equations  in 
order  to  describe  spacecraft  motion  about  a  reference  trajectory  for  the  purpose  of 
midcourse  guidance.  Before  solving  the  system,  he  transformed  to  cylindrical  coordinates 
and  to  chief  tme  anomaly  as  the  independent  variable  (except  for  the  out-of-plane 
component,  for  which  the  chiefs  eccentric  anomaly  became  the  independent  variable). 
Stem  solved  this  transformed  system  directly  in  closed  form.  He  then  derived 
expressions  for  the  in-plane  relative  position  coordinates  via  orbit-element  differences. 

He  also  showed  the  equivalence  of  the  two  methods,  relating  the  arbitrary  integration 
constants  in  the  original  solution  to  the  orbit-element  differences. 

The  next  contribution  to  the  third  path  was  from  Jones  in  1980  [46].  He  used  the 
time-domain  differential  equations  from  Stem,  and  presented  a  slightly  modified  form  of 
Stem's  closed-form  solution,  including  a  state  transition  matrix  (although  the  elements 
were  not  shown,  simply  the  algorithm  for  calculating  it).  He  also  demonstrated  that. 
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according  to  the  linearized  equations,  small  impulsive  changes  can  result  in  secular 
variations  (drift)  in  multiple  states,  not  only  along-track  position. 

Lawden  's  Equations:  Later  Developments. 

In  1987,  Carter  and  Humi  became  the  first  to  unify  two  of  the  paths  [14].  They 
recognized  that  Lawden's  equations  and  solution  are  equivalent  to  those  found  by  De 
Vries,  Tschauner  and  Hempel,  and  those  who  followed.  In  1990,  Carter  [II]  found  an 
improvement  to  the  definite  integral  in  Lawden's  solution.  Carter's  new  definite  integral 
is  no  longer  singular  for  circular  chief  orbits,  allowing  the  solution  to  be  applied  to  all 
chief  eccentricities. 

In  2002,  Inalhan,  Tillerson,  and  How  built  on  the  previous  results,  presenting 
time-domain  and  true-anomaly-domain  versions  of  Lawden's  equations  and  a  non¬ 
singular  version  of  the  solution  in  terms  of  six  integration  constants  [41].  They  showed 
that  setting  one  of  the  in-plane  integration  constants  to  zero  corresponds  to  a  boundedness 
condition;  that  is,  it  is  an  equivalent  constraint  to  requiring  the  chief  and  deputy  orbits  to 
have  the  same  semi-major  axis,  eliminating  the  secular  drift  terms.  One  of  their  most 
important  contributions  was  a  study  of  modeling  error  costs.  Specifically,  they  showed 
that  the  fuel  penalty  could  be  significant  for  using  an  HCW-based  formationkeeping 
control  algorithm,  even  for  Shuttle-type  eccentricities  (0.005). 

Also  in  2002,  Yamanaka  and  Ankersen  used  a  new  approach  to  eliminate  the 
definite  integral  in  the  solution  to  Lawden's  equations  [106].  Their  new  term  is 
proportional  to  the  transition  time;  thus  the  Yamanaka- Ankersen  solution  has  coefficients 
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expressed  in  terms  of  both  time  and  chief  true  anomaly.  They  organized  the  solution  into 
a  fairly  complex  state  transition  matrix. 

In  2003,  Broucke  solved  the  time-domain  version  of  Lawden's  equations  directly, 
without  having  to  transform  to  chief  true  anomaly  as  the  independent  variable  [8]. 

Instead,  he  replaced  the  time-varying  coefficients  first  found  by  De  Vries,  taking  their 
partial  derivatives  with  respect  to  the  chiefs  classical  orbit  elements.  He  showed  a 
method  for  organizing  the  solution  into  a  state  transition  matrix;  the  matrix  still  contains 
elements  that  vary  with  true  anomaly,  rather  than  explicitly  with  time. 

Later,  Sengupta  [83]  and  Sengupta  and  Vadali  [85,  86]  summarized  solutions  to 
the  elliptical  chief  problem.  They  solved  the  true-anomaly-domain  equations  of  motion 
in  terms  of  scaled  and  unsealed  coordinates,  using  the  Yamanaka-Ankersen  method  of 
eliminating  the  definite  integrals.  Their  "velocity"  solution  is  shown  first  in  the  true 
anomaly  domain  (i.e.,  they  solve  for  x'(t)  instead  of  i(t) ,  etc.).  They  included  a 

generally  applicable  mapping  from  the  solution's  six  constants  of  integration  to  Cartesian 
initial  conditions;  they  also  laid  out  an  approach  for  developing  a  state  transition  matrix  in 
these  terms.  Sengupta  and  Vadali  also  transformed  the  solution  into  a  time-explicit 
version,  retaining  the  infinite- series  Fourier-Bessel  expansions  [6]  of  the  coefficients; 
because  the  series  are  not  truncated,  they  remain  exact  solutions  to  Lawden's  equations, 
but  may  be  difficult  to  implement  practically.  Finally,  Sengupta  and  Vadali  developed 
simple  expressions  for  quantifying  the  amount  of  drift  in  unbounded-motion  cases  (i.e., 
mismatched  semi-major  axes);  they  showed  that  this  drift,  measured  over  an  integer 
number  of  chief  orbits,  is  a  minimum  when  that  measurement  occurs  at  chief  apogee. 
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In  2007,  Irvin  derived  a  new  form  of  Lawden's  equations  using  fractions  of  the 
chief  orbit  as  a  non-dimensional  independent  variable  [42] .  His  research  in  the  chief- 
orbit-fraction  domain  focused  on  applications  of  the  circular  chief  problem,  however. 

Equivalent  Forms. 

Other  linearized  relative  satellite  motion  models  for  the  elliptical  chief  problem 
are  physically  equivalent  to  using  Lawden's  equations,  even  though  they  are  not  described 
in  terms  of  Cartesian  LVLH  coordinates. 

In  1994,  Kelly,  aware  of  both  Lawden's  and  Tschauner  and  Hempel's  work, 
created  a  physically  equivalent  model,  but  resolved  in  the  coordinates  of  a  non-orthogonal 
reference  frame  [50].  In  Kelly's  frame,  the  2-direction  is  changed  so  that  it  is  always 
parallel  to  the  chiefs  inertial  velocity  vector.  Kelly  presented  time-domain  differential 
equations  of  motion  linearized  about  the  reference  orbit,  as  well  as  a  closed-form  solution 
in  the  form  of  a  relatively  simple  state  transition  matrix.  He  also  gave  the  mapping 
between  his  coordinates  and  more  traditional  orthogonal  LVLH  coordinates.  He  further 
showed  how  his  model  lends  itself  to  a  simple  impulsive  maneuver  scheme. 

In  2006,  Rathke  and  Izzo  used  a  reference  frame  defined  in  the  same  way  as 
Kelly's  to  describe  the  behavior  of  an  asteroid  after  an  impact  [72].  They  found  secular 
terms  only  in  the  2-direction.  They  asserted  that  their  solution,  when  rotated  into  LVLH 
coordinates,  is  "an  analytic  and  explicit  solution  of  the  Tschauner-Hempel  equations". 

Most  linearized  relative  satellite  motion  models  take  the  difference  between  the 
nonlinear  Keplerian  two-body  motion  of  each  satellite,  rotate  into  an  LVLH  frame,  and 
then  truncate  a  binomial  series  expansion  of  the  relative  gravitational  force  term.  In 
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1997,  Der  and  Danchick  took  a  different  approach  [25].  They  used  a  linearized  state 
transition  matrix  of  the  two-body  problem  itself  (using  universal  variables),  applying  it  to 
both  satellites  and  then  taking  the  difference.  This  resulted  in  a  relative  motion  solution 
in  ECI  coordinates,  which  can  be  transformed  into  other  reference  frames  via  rotation 
matrices. 

As  we  have  seen  from  Stern,  it  is  also  possible  to  represent  linearized  two-body 
motion  in  terms  of  small  differences  between  the  chief  and  deputy  orbit  elements.  This 
approach  is  physically  equivalent  to  using  Lawden's  equations,  but  yields  some  additional 
insight  into  how  the  large-scale  geometric  features  of  the  two  orbits  compare.  This 
concept  was  rediscovered  by  a  series  of  studies  that,  unaware  of  Stern's  results,  exploited 
orbit-element  differences  to  model  relative  satellite  motion. 

One  of  the  first  studies  to  do  so  was  that  of  Garrison,  Gardner,  and  Axelrad  in 
1995  [32].  They  derived  a  state  transition  matrix  and  LVLH  trajectory  equations 
linearized  with  respect  to  the  orbit-element  differences.  Their  numerical  results  were 
"nearly  identical"  to  a  numerical  integration  of  Lawden's  equations;  they  recognized  that 
any  differences  between  the  two  methods  would  have  to  be  of  second  order. 

The  best-known  work  using  orbit-element  differences  came  in  2003  and  2004 
from  Schaub  [77,  78]  and  Schaub  and  Junkins  [79].  This  research  used  a  set  of  elements 
in  which  eccentricity  and  argument  of  perigee  are  replaced,  so  that  the  system  is  not 
singular  for  circular  orbits.  Schaub  and  Junkins  included  a  linearized  mapping  between 
orbit-element  differences  and  Cartesian  coordinates,  and  they  derived  equations  for 
Cartesian  LVLH  velocity  by  directly  differentiating  the  position  equations.  They 
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examined  the  difference  in  mean  anomaly,  ,  as  a  function  of  time  and  of  chief  true 
anomaly,  and  derived  an  expression  for  drift  ( 5M  )  incorporating  the  effect  of  unequal 
orbit  energies  (mismatched  semi-major  axes). 

In  2006,  Lane  and  Axelrad  transformed  Broucke's  approach  to  develop  a  three- 
dimensional  linearized  description  of  relative  satellite  motion  in  terms  of  classical  orbit- 
element  differences  [54].  They  showed  that  the  in-plane  position  depends  on  the 
difference  in  right  ascension  of  the  ascending  node,  ,  in  a  way  that  Broucke  did  not 
capture,  because  Broucke  treated  coplanar  motion  independently  from  cross-track  motion. 

The  research  of  Sengupta  and  Vadali  [86]  also  included  orbit-element  differences. 
Using  classical  orbit  elements  as  well  as  a  variety  of  nonsingular  element  sets,  they 
derived  expressions  for  relative  position  and  velocity  in  terms  of  element  differences, 
showing  that  the  results  are  mathematically  equivalent  to  their  solutions  of  Lawden's 
equations.  They  included  the  mapping  from  the  six  constants  of  integration  to  these 
orbit-element  differences. 

Bounded-motion  Analysis  of  Lawden's  Equations. 

Some  of  the  most  interesting  and  geometrically  meaningful  analysis  of  the 
elliptical  chief  problem  has  occurred  in  studies  requiring  bounded  motion.  A  solution  to 
Lawden's  equations  can  be  made  bounded,  or  periodic,  if  it  is  subjected  to  one  of  the 
following  equivalent  constraints:  requiring  zero  drift,  equal  orbit  energies,  or  equal  semi¬ 
major  axes.  The  remaining  solution  terms  are  usually  somewhat  simpler. 

In  2003,  Zhang  and  Sun  offered  a  linearized  relative  satellite  motion  model  that 
assumes  equal  semi-major  axes  and  small  orbit-element  differences  [111].  They  gave 
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position  solutions  in  terms  of  these  differences  and  in  terms  of  Cartesian  initial 
conditions.  They  further  analyzed  the  geometry  of  the  resulting  trajectory  and  suggested 
formation  designs  not  possible  in  the  circular  chief  problem. 

Schaub's  research  [78,  79]  also  included  analysis  of  bounded  motion.  The 
linearized  expression  for  semi-major  axis  difference,  5a ,  can  be  set  to  zero  and  used  as  a 
boundedness  condition.  After  this  condition  is  imposed,  the  scaled  non-dimensional 
position  equations  of  Schaub's  model  give  some  geometrical  insight  into  the  trajectory. 

The  research  of  Lane  and  Axelrad  [54]  resulted  in  a  simple  bounded-motion 
model  with  parameters  very  similar  to  the  ROEs  of  the  circular  chief  problem.  The  eight 
parameters  are  not  linearly  independent,  and  they  are  defined  in  terms  of  five  orbit- 
element  differences  (only  five  because  5a  is  already  set  to  zero  in  the  boundedness 
condition).  They  analyzed  three  specific  formation  geometries  possible  in  the  bounded 
motion  model.  The  most  interesting  is  an  in-track/cross-track  formation,  in  which  the 
deputy's  motion  follows  an  ellipse  in  the  y-z  plane  with  no  associated  x-motion,  behavior 
impossible  in  the  circular  chief  problem. 

The  research  of  Sengupta  and  Vadali  [86]  also  included  geometrically  insightful 
bounded-motion  analysis.  They  derived  a  boundedness  condition  from  their  solution  to 
Lawden's  equations  by  setting  one  of  the  constants  of  integration  to  zero;  they  expressed 
this  condition  in  terms  of  both  scaled  and  unsealed  Cartesian  coordinates.  They  next 
wrote  a  bounded  version  of  the  position  solution  using  five  design  parameters;  again, 
there  are  only  five  design  parameters  because  one  degree  of  freedom  was  used  to 
eliminate  drift.  After  analyzing  the  effect  of  these  design  parameters  and  comparing  to 
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the  bounded  HCW  model,  Sengupta  and  Vadali  observed  five  geometric  effects  of 
eccentricity:  higher-frequency  harmonics,  amplitude  scaling,  phase  shift,  off-center 
formations,  and  skewness  of  the  relative  orbit  "plane." 

In  2008,  Jiang,  Li,  Baoyin,  and  Gao  used  a  unique  celestial-sphere  approach  to 
derive  a  relative  satellite-motion  model  [45].  They  showed  that  in  its  linearized  and 
bounded-motion  form,  their  model  is  equivalent  to  a  bounded-motion  solution  to 
Lawden's  equations.  They  created  a  position  solution  using  five  parameters,  relating 
these  parameters  to  the  integration  constants  used  by  Inalhan,  Tillerson,  and  How,  and 
then  to  Cartesian  initial  conditions.  Next,  they  used  a  change  of  variable  to  ^(t) ,  defined 
as 


'  '  2 

In  this  way,  they  eliminated  all  trigonometric  terms,  creating  a  completely 
algebraic  parameterized  solution.  This  algebraic  form  lends  itself  to  simple  geometric 
analysis;  Jiang  et  al.  found  that  the  trajectory's  projection  into  each  of  the  LVLH 
coordinate  planes  could  be  described  by  a  quartic  equation.  They  also  found  that  such  a 
bounded  trajectory  could  have  at  most  one  self-intersection  point;  they  derived 
expressions  for  locating  such  points,  as  well  as  for  locating  chief-deputy  collision  points 
(i.e.,  when  the  deputy  crosses  the  LVLH  frame's  origin).  They  further  proved  that,  except 
in  a  limited  number  of  degenerate  cases,  a  bounded  relative  trajectory  is  not  confined  to  a 
plane,  as  is  the  case  in  the  circular  chief  problem. 
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The  most  dramatic  result  from  Jiang  et  al.  was  proving  that,  in  three  dimensions,  a 
bounded  relative  trajectory  is  always  confined  to  a  quadric  surface.  The  coefficients  of 
the  quadric-surface  equation  can  be  related  to  chief  eccentricity  and  the  five  relative  orbit 
parameters.  Jiang  et  al.  derived  expressions  for  finding  the  geometric  center  of  the 
quadric  surface;  for  rotating  from  LVLH  coordinates  to  quadric-surface  principal- axis 
coordinates;  and  for  determining  in  a  particular  case  whether  the  quadric  surface  will  be 
an  elliptic  cylinder,  a  hyperboloid  of  one  sheet,  or  an  elliptic  cone.  All  self-intersecting 
relative  trajectories  are  confined  to  the  surface  of  an  elliptic  cone. 

Low-eccentricity  Approximations. 

Since  relative  satellite  motion  models  for  the  elliptical  chief  problem  tend  to  be 
somewhat  complex,  there  have  been  numerous  efforts  to  simplify  the  equations  of  motion 
or  their  solutions  by  assuming  small  values  for  the  chiefs  eccentricity.  Many  of  these 
expand  coefficients  containing  the  chiefs  true  anomaly,  angular  rate,  and  orbit  radius  in 
an  eccentricity  power  series,  truncating  the  results  to  retain  only  terms  which  are  first-, 
second-,  or  third-order  in  eccentricity.  The  more  terms  retained,  the  more  complex  the 
result,  but  also  the  more  accurate  for  a  given  level  of  chief  eccentricity.  This  series- 
expansion  approach  can  have  the  added  operational  advantage  of  rendering  the  solution 
completely  explicit  in  time. 

In  1965,  Anthony  and  Sasaki  solved  their  second-order  equations  of  relative 
motion  via  the  method  of  differential  corrections  [4].  That  is,  their  relative  motion 
solution  was  composed  of  the  HCW  solution  plus  a  correction  for  first-order  eccentricity 
plus  another  correction  for  quadratic  terms  in  the  displacement.  Leaving  out  the 
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eccentricity  correction  resulted  in  a  solution  equivalent  to  London's.  But  leaving  out  the 
quadratic  correction  resulted  in  a  linear  model  valid  for  first-order  eccentricity.  This 
linear  model  is  time-explicit,  containing  periodic  terms  with  frequencies  and  2n^t 

(where  is  the  chiefs  mean  motion)  and  mixed-secular  terms.  In  numerical  examples, 

the  linear  model  did  not  perform  as  well  as  the  quadratic,  as  expected,  especially  when  it 
came  to  capturing  drift  in  the  radial  direction;  but  the  example  formations  were  all  on  the 
scale  of  50  miles  or  larger. 

Jones  [46]  also  addressed  small-eccentricity  orbits  in  his  research,  using  a  first- 
order  truncated  eccentricity  power  series.  He  gave  the  elements  of  the  resulting  state- 
transition  matrix,  and  showed  that  the  percentage  error  using  this  approximation  is 
significantly  improved  over  that  using  a  circular  orbit  approximation  (i.e.,  neglecting  the 
chief  eccentricity  and  using  HCW  dynamics  alone). 

In  2000,  Melton  found  an  approximate  solution  to  the  time-domain  version  of 
Lawden's  equations  in  both  rectangular  and  cylindrical  coordinates  [66].  Using 
Lagrange's  generalized  expansion  theorem  [6],  he  expanded  the  dynamics  matrix,  creating 
a  time-explicit  matrix  equation  of  motion  in  the  form  of  an  eccentricity  power  series 
which  can  be  truncated  at  the  desired  order.  This  form  admits  a  solution  in  the  form  of  a 
state  transition  matrix  expanded  in  an  eccentricity  power  series,  truncated  to  the  same 
order.  The  zeroth-order  matrix  corresponds  to  HCW  dynamics.  Melton  listed  the 
elements  for  the  first-  and  second-order  state  transition  matrices  in  both  rectangular  and 
cylindrical  coordinates.  Of  course,  the  method  described  could  be  used  to  generate 
approximate  solutions  of  arbitrarily  high  order.  The  first-  and  second-order  elements 
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include  periodic  terms  with  frequencies  ,  2n^t ,  and  ,  plus  mixed- secular  terms. 


Melton  believed  that  the  mixed-secular  terms  in  the  x-direction  lacked  physical 
significance. 

In  2003,  Sabol,  McLaughlin,  and  Luu  proposed  a  new  relative  satellite  motion 
model,  which  they  entitled  the  COWPOKE  Equations,  standing  for  Cluster  Orbits  with 
Perturbations  of  Keplerian  Elements  [76].  In  spite  of  the  title,  the  original  model  did  not 
include  perturbing  forces,  making  it  an  approximation  of  the  two-body  elliptical  chief 
problem.  The  model  used  a  linearized  orbit-element  difference  approach  to  describe 
relative  position  in  spherical  coordinates;  it  then  used  Eourier-Bessel  series  expansions  to 
make  the  equations  explicit  in  time,  allowing  truncation  at  the  desired  power  of  chief 
eccentricity.  They  carefully  analyzed  the  question  of  how  many  terms  to  retain  in  order  to 
approximate  well  even  for  higher  eccentricities;  this  is  important  because  for  values 
higher  than  about  0.03,  the  coefficients  may  increase  faster  than  the  powers  of 
eccentricity  decrease.  Sabol  et  al.  wrote  out  some  of  the  required  expansions  to  tenth 
order  and  the  complete  equations  to  first  order.  Their  numerical  simulations  showed 
reasonable  results  for  the  first-order  model  at  a  chief  eccentricity  of  0.01  and  for  the 
tenth-order  model  at  a  chief  eccentricity  of  0.7. 

In  2003,  Vaddi,  Vadali,  and  Alfriend  extended  previous  approximations  by  adding 
corrections  to  the  HCW  solutions  that  are  quadratic  in  the  relative  distance  and  third- 
order  in  the  chiefs  eccentricity  [97].  A  linearized  model  is  easily  constructed  by  leaving 
out  the  quadratic  corrections;  Vaddi,  et  al.,  wrote  out  the  linear  equations  in  terms  of 
Cartesian  initial  conditions  for  the  special  case  when  motion  is  bounded  and  epoch  time 
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occurs  at  chief  perigee.  Interestingly,  the  out-of-plane  equation  contains  no  third-order 
terms.  They  also  included  a  boundedness  condition  for  this  special  case. 

Schaub's  research  [78,  79]  also  examined  a  small-eccentricity  approximation. 
Noting  that  the  linearization  of  his  orbit-element  difference  method  already  assumed  a 
small  value  for  the  ratio  of  relative  distance  to  chief  orbit  radius,  he  stated  that  powers  of 
chief  eccentricity  smaller  than  that  ratio  could  be  neglected.  He  created  a  first-order 
approximation  of  chief  orbit  radius  (see  Eq.  (D.2))  as  (l  -  cos  (t)) ,  where  the 

subscript  C  indicates  chief  orbit  elements.  He  was  then  able  to  write  first-order 
approximations  for  SM  and  for  the  relative  position  coordinates. 

In  2006,  Ketema  developed  a  complex  nonlinear  relative  satellite  motion  model, 
and  then  created  a  linearized  approximation  by  assuming  small  chief  eccentricity  and 
small  differences  in  eccentricity  and  inclination  [51].  After  further  assuming  bounded 
motion,  Ketema  created  a  model  with  interesting  geometric  characteristics.  In  this  model, 
the  relative  trajectory's  projection  into  the  chief  orbit  plane  (the  x-y  plane)  is  shown  to  be 
an  ellipse  centered  at  (o,  -a^,  +  Se)  sin  o )  ’  where  ^  is  the  chiefs  eccentric 

anomaly  at  epoch. 

Quadratic  Models. 

As  with  the  circular  chief  problem,  a  relative  satellite  motion  model  that  includes 
quadratic  terms  in  the  relative  distance  increases  accuracy  over  a  linearized  model, 
especially  as  the  relative  distance  grows. 

The  model  of  Anthony  and  Sasaki  [4]  was  the  first  such  model.  Their  quadratic 
corrections  include  periodic  terms  with  frequencies  of  n^t  and  2n^t ,  mixed-secular 
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terms,  and,  in  the  along-track  direction  only,  terms  with  a  factor  of  .  They  also  showed 
how  this  model,  continuing  with  the  method  of  differential  corrections,  can  be  used 
effectively  to  solve  rendezvous  problems. 

In  1967,  Euler  and  Shulman  derived  true-anomaly-domain  equations  of  relative 
motion  directly  analogous  to  Lawden's  equations,  but  retaining  quadratic  gravity  terms 
[28].  To  solve,  they  used  the  known  solutions  to  Lawden's  equations  as  the  homogeneous 
solution,  converted  the  quadratic  terms  to  forcing  functions,  and  solved  for  the  particular 
solutions  in  the  x-  and  z-directions  via  variation  of  parameters.  They  were  unable  to  find 
a  closed-form  particular  solution  in  the  y-direction. 

The  model  developed  by  Vaddi,  Vadali,  and  Alfriend  [97]  included  quadratic 
corrections,  as  we  have  seen.  They  also  analyzed  the  higher-order  coupling  between 
nonlinearity  and  eccentricity,  although  this  coupling  was  not  reflected  in  their  model. 

In  2006,  Sengupta,  Sharma,  and  Vadali  extended  the  results  of  Euler  and  Shulman 
by  finding  a  closed- form  solution  in  all  three  relative  components  [84].  They  were  able  to 
do  so  in  the  case  of  bounded  motion.  To  solve,  they  used  a  perturbation  approach 
involving  a  small  parameter  which  depends  on  relative  distance  and  chief  eccentricity. 
Since  bounded  motion  was  important  to  their  result,  they  derived  an  index  to  quantify 
drift,  or  rather,  to  quantify  error  as  compared  with  the  boundedness  condition.  They  also 
noted  that,  strictly  speaking,  all  two-body  solutions  are  bounded,  but  that  what  we  really 
need  is  1:1  resonance  between  the  chief  and  deputy  orbit  periods. 
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Applications  of  Relative  Motion 

Relative  satellite  motion  models  have  been  used  for  a  variety  of  spaceflight 
applications.  In  1981,  Kelley,  Cliff,  and  Lutze  used  HCW  dynamics  and  a  maneuver- 
scheme  game  scenario  to  analyze  pursuit  and  evasion  on  orbit  [49] .  As  we  have  seen, 
Rathke  and  Izzo  [72]  applied  a  linearized  elliptical  chief  model  to  study  asteroid 
deflection  missions. 

Some  studies  have  applied  a  variety  of  models  to  develop  relative  guidance, 
navigation,  and  control  applications  [22,  27,  34,  47,  60,  63,  87,  89,  96,  100,  109]. 

Other  studies  have  examined  formation  flying  missions  [9,  35,  63,  75,  95,  107, 
111],  such  as  the  European  PRISMA  mission  [22],  NASA's  Magnetosphere  Multiscale 
mission  [21,  34,  62,  80,  109],  or  a  variety  of  Earth-observing  synthetic-aperture  or  bistatic 
radar  concepts  [30,  52,  81]. 

The  largest  group  of  studies  considered  for  this  survey  considered  rendezvous  and 
proximity  operations  (RPO)  [7,  15,  31,  38,  42,  43,  44,  69,  83,  92,  104,  105,  110]. 

Specific  studies  have  covered  operations  near  the  International  Space  Station,  including 
Eree  Elyer  missions  [64,  74]  and  the  European  ATV  resupply  vehicle  [89,  99].  Early 
autonomous  RPO  experiments  have  included  the  Japanese  ETS-Vn  mission  [68],  the  Air 
Force’s  XSS-11  mission  [104],  and  Orbital  Express,  managed  by  DARPA  and  NASA 
[29,  100].  In  2005,  Naasz  planned  for  a  robotic  servicing  mission  to  the  Hubble  Space 
Telescope  by  studying  a  parameterized  version  of  the  HCW  model,  using  a  variety  of 
geometric  constraints  [67] . 
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Higher-fidelity  Models 

The  present  research  focuses  on  linear  or  quadratic  approximations  of  unperturbed 
two-body  relative  satellite  motion.  However,  for  a  model  to  be  accurate  over  long 
distances  or  long  time  scales,  it  must  consider  a  variety  of  other  factors.  The  literature 
contains  many  examples  of  such  high-fidelity  models. 

High-order  nonlinear  models  can  be  accurate  even  without  a  closeness  assumption 
that  requires  the  deputy  to  remain  within  a  certain  distance  from  the  chief  [17,  18,  19,  20, 
37,  47,  51,  53,  70,  87,  111].  These  models  tend  to  be  complex,  and  several  require 
determining  both  the  chief  and  deputy  orbit  elements  prior  to  calculating  the  relative 
motion. 

In  2009,  Segal  and  Gurfil  developed  a  nonlinear  six-degree-of-freedom  model  that 
included  spacecraft  rotation  effects  [82]. 

Many  models  consider  the  differential  effects  of  the  J2  [3,  23,  33,  39,  50,  77,  78, 
81,  83,  108]  term  of  the  Earth's  gravity  field.  Wiesel  accounted  for  higher-order  gravity 
terms,  as  well  [103].  Other  models  considered  the  differential  effects  of  air  drag  in  low- 
altitude  orbits  [5,  13,  23,  77,  103].  Lovell,  Horneman,  Tragesser,  and  Tollefson 
incorporated  the  J2  and  differential  drag  effects  on  ROEs  [59]. 

Einally,  even  though  the  present  research  does  not  consider  eccentricity  greater 
than  or  equal  to  1,  some  researchers  have  considered  relative  satellite  motion  on 
parabolic,  hyperbolic,  or  rectilinear  Keplerian  orbits  [10]. 


31 


Model  Comparisons 

A  few  studies  in  the  literature  have  been  devoted  to  comparing  relative  satellite 
motion  models  or  to  developing  criteria  for  such  comparisons. 

For  example,  in  1998,  Carter  provided  a  fairly  comprehensive  bibliography  and 
taxonomy  of  models  developed  to  that  point  [12].  He  was  one  of  the  few  authors  aware 
of  the  third  path  of  development  in  the  first  generation  of  linearized  elliptical  chief 
models  (Stern  [90]  and  Jones  [46]). 

In  2003,  Melton  numerically  compared  four  models  to  HCW  and  to  numerically 
integrated  fully  nonlinear  equations  of  motion  [65].  The  models  he  compared  are 
Yamanaka-Ankersen  [106],  Broucke  [8],  Melton  [66],  and  Karlgaard-Lutze  [48].  The 
first  two  performed  the  same;  this  is  understandable,  since  both  are  versions  of  Lawden's 
equations.  Except  for  chief  eccentricities  below  about  0.2,  the  first  two  performed 
noticeably  better  than  the  third;  this,  too,  is  understandable,  since  Melton's  model  is  an 
approximation  (second-order,  in  this  case)  to  Lawden's  equations.  The  fourth  model 
performed  much  better  than  HCW,  also  as  expected. 

In  2005,  Alfriend  and  Yan  showed  that  the  accuracy  of  a  relative  satellite  motion 
model  increases  with  its  complexity  [2] .  They  compared  a  variety  of  models  and 
discovered  that  the  scale  of  the  relative  trajectory  and  the  chief  eccentricity  are  key 
parameters  for  characterizing  their  accuracy.  Their  evaluation  considered  not  only 
eccentricity,  but  also  J2  and  high-order  nonlinear  effects. 

Also  in  2005,  Lovell,  Chavez,  and  Johnson  compared  the  Yamanaka-Ankersen 
model  to  HCW  dynamics  and  to  numerically  simulated  fully  nonlinear  two-body  motion 
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[58].  They  found  that  for  the  examples  considered,  for  chief  eccentricity  about  0.001,  the 
HCW  model  was  significantly  less  accurate.  They  also  found  that  the  accuracy  of  the 
Yamanaka-Ankersen  model  began  to  drop  off  at  values  of  chief  eccentricity  above  about 
0.6.  Since  any  version  of  Lawden's  equations,  including  Yamanaka-Ankersen,  is  a 
linearized  model,  it  cannot  capture  the  high-order  coupling  between  eccentricity  and  the 
nonlinear  terms  [97];  it  is  possible  that  this  effect  could  explain  the  results  from  Lovell  et 
al. 

Summary 

One  goal  of  this  thesis  is  to  analyze  the  ROE  realization  of  the  HCW  model  for 
the  circular  chief  problem,  identifying  some  of  the  features  which  provide  operational 
efficacy.  The  other  main  goal  is  to  find  a  new  model  for  the  elliptical  chief  problem 
which  retains  the  advantages  of  ROEs. 

After  searching  the  extensive  literature  for  such  a  model,  we  found  that  none  had 
all  of  the  following  characteristics:  simple  parameterized  trajectory  equations,  parameters 
with  geometrically  intuitive  interpretations,  validity  for  eccentric  chief  orbits,  and  validity 
for  drifting  relative  trajectories. 

Many  existing  models  exhibit  one  or  more  of  these  characteristics,  but  none 
exhibit  all.  Chapter  IV  will  explore  several  possible  approaches  for  developing  a 
satisfactory  model. 
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III.  Relative  Orbit  Elements  for  the  Circular  Chief  Problem 


As  we  have  seen,  the  most  commonly  applied  relative  satellite  motion  model  is 
the  Hill-Clohessy-Wiltshire  (HCW)  model  for  the  circular  chief  problem,  and  Relative 
Orbit  Elements  (ROEs)  are  a  particularly  useful  realization  of  that  model.  The  goal  of 
this  chapter  is  to  describe  the  characteristics  of  ROEs  which  provide  such  a  high  degree 
of  geometrical  insight  and  operational  efficacy.  After  an  intuitive  derivation  of  the  ROE 
parameters,  we  will  also  see  that  the  instantaneous  trajectory  of  the  deputy  satellite  is  an 
ellipse,  viewed  in  three-dimensional  space. 

The  HCW  Model 

Derivations  and  solutions  of  the  HCW  equations  of  motion  can  be  found  in  many 
astrodynamics  texts  [79,  98,  102].  Assuming  no  external  forces  (except,  of  course, 
inverse-square  gravity),  the  scalar  equations  of  motion  are 


X  -  2n^y  -  3njx  =  0 

y  +  2n^x  =  0  (3.1) 

Z  +  nj'z  =  0 

where  is  the  chiefs  mean  motion.  In  the  circular  chief  problem,  is  also  equal  to  its 
angular  rate  at  every  point  in  the  orbit.  Using  the  notation  introduced  in  Chapter  I,  the 
angular  velocity  of  the  chiefs  LVLH  frame  with  respect  to  inertial  space  is 
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=  0 


Note  that  the  z-motion  is  uncoupled  from  the  in-plane  motion. 

If  we  use  the  relative  position  and  velocity  at  epoch  time  as  initial  conditions, 
denoted  by  the  subscript  0,  then  the  position  solution  to  the  HCW  equations  is 


x(t)  =  —  sinTn^  -(3vq  +  ^^)cosrnp  +  4vg 

4  *  2jc 

y(t)  =  (6vq  h - ^)sin[n(,  -l - ^cos[np  (t-to)]-(6n^VQ  -f3>’o)(t-tg) 

/T  /-» 


lit)  =  Zg  cos  [n^  (t  -  to )]  +  —  sin  [«(,(?-  to )] 


Again  using  the  notation  of  Chapter  I  for  relative  position  and  velocity,  the 
position  can  be  expressed  in  terms  of  state  transition  sub-matrices,  so  that 


(Cto)^Z^  (to)  +  ®„  (^’^o) 


where 
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-3cos[«^(f-fQ)]  +  4  0  0 


sin[?i^(f-f„)] 


0  cos[«^(f-f„)] 

2cos[n^{t-t^)]  ^  2 


2cos[n^{t-t^)]  2  4sin[«^(f-fJ] 


-3(t-tJ  0 


sin[«^(f-fo)] 


(3.4) 


Likewise,  the  velocity  solution  can  be  written  as 


d  ^  d 

^  ^  (0  “  ®  vr  (^’  ^0  )  ^  (^0  )  +  ®  w  (^’  ^0  )  ^  ^  (^0  ) 


(3.5) 


where 


sin[n(^  (t- 

»«)]  0 

0 

®v.(ri^o)  = 

6^^  cos[np  {t-tfj) 

]-6np  0 

0 

0 

0  -n^  sin 

[nc{t- 

^o)]_ 

cos[nc(t-to)] 

2sin[nc(t-to)] 

0 

e 

o 

II 

4cos[n^{t-t,)]- 

3 

0 

0 

0 

cosj 

nc{, 

(3.6) 


A  Geometric  Parameterization  of  the  Solution 

This  section  is  a  detailed  derivation  of  the  six  ROEs:  (t) ,  ,  (5^  or  y0(t) , 

Zmax  ’  Vq  or  Y  [60].  To  characterize  the  geometry  of  the  HCW  relative  trajectory,  we 
must  examine  the  various  terms  in  Eqs.  (3.2).  The  x-equation  is  periodic,  with  the  third 
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To  deal  with  the  periodic  terms  in  the  x-,  y-,  and  z-equations,  we  must  use  the 
following  principle,  the  Harmonic  Addition  Theorem.  (See  Appendix  B.)  A  periodic 
function  defined  as 
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/(^)  =  Asin^  +  ficos^ 


(3.9) 


can  be  equivalently  written  both  as 

f(d)  =  4^7Wsm(d  +  atan2(B,A))  (3.10) 


and  as 


f(0)  =  -^Ia^  +B^  cos(0  +  atan2(A,-B)) 


(3.11) 


where  atan2  is  the  quadrant- specific  inverse  tangent  function.  Both  of  these  forms,  Eqs. 

(3.10)  and  (3.11),  are  periodic  functions  whose  amplitude  (  +B^  ),  initial  phase 

(atan2{B,A)  or  atan2{A,-B) ),  and  time-varying  phase  parameter  (^ )  are  all  readily 
apparent. 

The  periodic  terms  of  the  y-equation  from  Eqs.  (3.2)  can  be  re-written  using  Eq. 

(3.10) ,  yielding 


y(t)  = , 


6xq  + 


4^0 


n, 


+ 


c  J 


f 

J 


sm 


{t  -  tff)  +  atan2 


V 


n, 


c  J 


+  yd 


The  amplitude  of  the  y-oscillation  can  be  written  as  ,  another  ROE: 


a,  =  2. 1  3vn  -I- 


2jo 


n. 


c  J 


f  .  3^ 


V^c  y 
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Also,  factoring  —  out  of  both  arguments  of  the  atan2  function,  we  can  describe  a 

ric 

constant  phase  parameter 


Pq  =  +  2yP 


We  can  use  P^  as  one  of  the  ROEs,  or  instead  we  can  define  a  time-varying  phase 
parameter 


/?  =  y?o+'^c(^-^o) 


(3.12) 


Then  we  can  describe  the  along-track  motion  as 


y(t)  =  a^smp  +  y^ 


(3.13) 


Next,  we  can  use  Eq.  (3.11)  to  re-write  the  x-equation  from  Eqs.  (3.2),  yielding 


.^(0  = 


\^c  J 


-h 


3X(,  -I- 


2jo 


cos 


‘c  J 


He  {t  -tP)  +  atanl 


^,3^0  + 

V^C 


2jo 


‘c  J 


+  x. 


Note  that  the  amplitude  of  the  x-oscillation  equals  .  Also,  multiplying  both 

arguments  of  the  atan2  function  by  ,  we  see  that  the  x-motion  is  90  degrees  out  of 
phase  with  the  y-motion,  such  that 
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a  ^ 

x(t)  =  -^cos  p  +  x^ 


(3.14) 


Finally,  applying  Eq.  (3.10)  to  Eq.  (3.2),  we  can  re-write  the  z-equation  as 


(7  ^ 

2 

- 

(  •  3 

sin 

He  {t  -1^)  +  atanl 

z  — 

^0, 

v^c  y 

V  y 

Thus,  the  deputy's  cross-track  motion  is  a  simple  sinusoid  motion  with  amplitude  , 

first  identified  as  a  useful  parameter  by  Clohessy  and  Wiltshire  [16],  and  adopted  later  as 
another  ROE: 


z 


max 


-I- 


r  •  3 


v^c  y 


(3.15) 


Eor  the  sixth  and  final  ROE,  we  could  define  a  constant  phase  parameter 


=  atanl 


( 


‘c  ) 


(3.16) 


Or,  alternatively,  we  could  use  the  constant  phase  difference  between  the  cross-track 
oscillation  and  the  along-track  oscillation: 

y  =  ¥o-Pi, 


Thus,  the  z-motion  is 
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=  siniP  +  y) 


(3.17) 


We  have  now  used  ROEs  to  describe  the  deputy's  relative  trajectory  in  three 
parametric  equations  (Eqs.  (3.13),  (3.14),  and  (3.17)),  as  did  Eovell  and  Tragesser  [60]. 

The  Three-dimensional  Geometry  of  the  Relative  Trajectory 
Projection  into  the  Orbit  Plane. 

To  understand  the  physical  geometry  of  the  deputy's  motion,  as  seen  from  the 
chief,  we  look  first  at  the  trajectory's  projection  into  the  chiefs  orbit  plane,  the  x-y  plane. 
Eirst  let  us  consider  the  degenerate  case  when  equals  zero.  The  deputy's 

relative  position  is  simply  (x^,  (t)) ,  meaning  the  projected  trajectory  is  simply  a  point, 

possibly  drifting  in  the  y-direction. 

If  is  not  zero,  we  can  re-write  Eq.  (3.13)  as 

y(0-y^  =sin/?  (3.18) 

ae 


and  then  square  both  sides,  obtaining 


2 


sin^  P 


(3.19) 


Next,  rewriting  Eq.  (3.14)  as 
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(3.20) 


x{t)-x^ 

/2 


-cos  P 


and  again  squaring  both  sides  yields 


{x{t)-xPf 


=  COS^  P 


(3.21) 


Adding  Eqs.  (3.19)  and  (3.21)  and  then  applying  the  trigonometric  identity  in  Eq. 

(A.l), 


(x(0-xj^ 


+ 


„  2 


=  1 


(3.22) 


This  is  the  equation  of  an  ellipse  in  the  x-y  plane  centered  at  point  ( ),  a 
point  which  is  drifting  because  of  the  secular  term  in  the  definition  of  ,  Eq.  (3.8).  The 


ellipse's  semi-major  axis  is  and  lies  along  the  y-axis;  the  semi-minor  axis  is 


and 


lies  along  the  x-axis.  Thus  the  deputy's  in -plane  motion  can  be  described  as  an 
instantaneous  2x1  ellipse. 

Instantaneous  Plane  of  the  Trajectory. 

To  understand  the  trajectory  in  three  dimensions,  I  will  follow  an  argument  used 
by  Press  [71]  to  describe  stationary  satellite  formations  (in  which  the  deputy's  along-track 
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secular  drift  is  constrained  to  be  zero).  Here,  I  will  extend  Press's  result  by  describing  the 
general  case,  with  drift  allowed. 

We  have  already  used  Eq.  (3.22)  to  describe  an  instantaneous  ellipse  in  the  x-y 
plane.  Viewed  in  three  dimensions,  however,  Eq.  (3.22)  describes  an  elliptical  cylinder 
orthogonal  to  the  x-y  plane.  Thus,  at  any  given  time,  the  deputy's  relative  trajectory  is 
confined  to  the  surface  of  this  elliptical  cylinder. 

Next,  we  can  use  the  trigonometric  identity  in  Eq.  (A. 2)  to  expand  Eq.  (3.17), 
obtaining 


z{t)  =  (cos  p?,my  +  sin  cos  y) 


(3.23) 


We  can  then  substitute  Eqs.  (3.18)  and  (3.20)  into  Eq.  (3.23)  to  obtain 


z(t)  =  z^ 


x{t)-Xj.  .  y{t)-y. 

, mi y  +  — —cosy 


Rearranging, 


0  =  (-2  sin  y){xit)  -  x J  +  (cos  y)(yit)  -y^)  — —  {zit)  -  0)  (3 .24) 

7 

^max 

This  is  the  equation  of  a  plane  passing  through  the  point  (x^,y^  ,0),  the  center  point  of  the 
x-y  projected  ellipse.  The  plane  has  normal  vector 
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(3.25) 


Thus  at  any  time,  Eq.  (3.24)  must  be  true,  and  the  deputy's  motion  must  be 
confined  to  the  plane  it  describes.  This  plane  is  constant  in  orientation,  since  is  a 
constant  vector,  but  it  drifts  along  with  the  center  point  of  the  x-y  projected  ellipse.  Note 
that  the  z-component  of  N  is  always  negative;  thus  N  always  points  “below”  the  chiefs 
orbit  plane. 

The  deputy's  relative  trajectory  in  3-D  must  therefore  be  on  the  intersection  of  this 
plane  with  the  elliptical  cylinder  described  by  Eq.  (3.22).  The  orientation  of  this  plane 
with  the  c-frame  can  be  easily  understood  from  Eq.  (3.25).  In  fact,  the  3-D  instantaneous 
trajectory  and  its  drift  can  be  completely  described  by  the  chiefs  mean  motion  (n^)  and 
the  five  non-periodic  ROEs  of  the  deputy's  relative  motion  (constants  ,  and 

y ,  and  secular  element  ).  The  deputy's  motion  within  its  3-D  trajectory  can  be 
completely  described  by  the  sixth  ROE,  p . 

The  True  Ellipse. 

We  know  (see  Appendix  C)  that  any  trajectory  equations  which  satisfy  the 
following  form  must  describe  an  ellipse: 
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x{t)  =  au^  cos  6  +  bv^  sin  0  +  C^ 
y(t)  =  au,,  cos 6  +  bv^,  sin 6  +  C, 
z(t)  =  au^  cos  6  +  bv^  sin  0  +  C^ 


(3.26) 


If  we  can  show  that  this  form  is  satisfied  by  the  deputy's  c-frame  parametric 
equations,  then  we  will  have  proven  that  the  deputy's  relative  trajectory  is  an 
instantaneous  ellipse.  However,  we  must  consider  carefully  what  form  of  the  parametric 
equations  to  use. 

The  derivation  in  Appendix  C  uses  a  and  b  as  the  semi-major  and  semi-minor 
axis  lengths,  and  it  defines  the  first  two  basis  vectors  of  the  principal-axis  reference  frame 
as  aligned  with  the  ellipse's  major  and  minor  axes.  However,  it  makes  no  statement  about 
whether  the  major  axis  corresponds  to  a  and  the  reference  frame's  1 -direction  or  to  b 
and  the  reference  frame's  2-direction.  All  we  know  is  that  is  a  phase  parameter  moving 
in  a  positive  rotation  from  the  1 -direction  to  the  2-direction,  as  shown  in  Figure  2. 


Figure  2.  Direction  of  theta 
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If  we  wish  to  map  correctly  into  the  ellipse's  principal- axis  frame  from  the  c- 
frame,  we  must  use  an  appropriate  phase  parameter.  /?  is  not  ideal  because  it  rotates  in  a 
negative  sense  starting  from  the  negative  x-axis,  as  shown  in  Figure  3.  To  prove  this,  let 
us  introduce  a  new  definition  of  fi ,  equivalent  to  Eq.  (3.12).  Replacing  the  initial 
conditions  within  the  definition  of  [3^  with  time-varying  components,  and  assuming  for 
illustration  purposes  that  is  zero,  we  see  that 

/?  =  atan2(x,  3n^x  -f  2y) 


(along- 


Figure  3.  Direction  of  beta 
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We  know  that  [3  is  rotating  negatively  (clockwise  in  Figure  3)  because  when  the 


deputy  is  at  perigee,  it  is  at  a  lower  altitude  than  the  chief.  Thus,  it  is  moving  faster  than 
the  chief,  and  the  y-coordinate  of  its  position  must  be  increasing;  in  other  words,  y  must 
be  greater  than  zero.  The  opposite  is  true  at  apogee. 

Where  we  begin  measuring  the  rotation  of  ,  it  must  take  on  a  zero  value.  This 
occurs  where  the  sine  of  x  equals  zero;  in  other  words,  where  x  equals  zero.  This  is  true 
at  the  maximum  and  minimum  values  of  x;  examining  the  2x1  projected  ellipse  in  the  x-y 
plane  makes  it  clear  that  these  values  occur  at  apogee  and  perigee,  respectively.  To 
determine  where  [5  begins  its  rotation,  we  must  find  a  way  to  discriminate  between  these 
two  possible  points. 

A  phase  angle  moves  through  its  first  quadrant  where  both  arguments  of  its  atan2 
function  are  positive,  just  after  the  phase  angle  equals  zero.  moves  through  its  first 
quadrant  where  both  the  sine  of  x  and  the  cosine  of  3n^x  +  2y  take  on  positive  values; 

we  must  determine  whether  this  is  true  just  after  perigee  or  just  after  apogee.  The  sine  of 
X  near  zero  is  positive  only  where  x  is  positive;  that  is,  where  x  is  increasing.  This  is 
not  true  in  the  quadrant  clockwise  from  the  apogee  point;  it  is  true  in  the  quadrant 
clockwise  from  the  perigee  point,  as  shown  in  Figure  3. 

Thus,  P  rotates  in  a  negative  sense,  starting  from  the  negative  x-axis. 

To  address  the  concern  about  the  direction  of  rotation,  we  can  constrain  the  third 
principal  axis  w  to  have  an  opposite  sense  from  the  third  LVLH  axis,  ,  so  that 
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N 

w  =  — — 

N 

2 

Still,  is  a  phase  variable  measured  from  the  negative  x-axis,  while  0  from  Eqs. 

(3.26)  is  a  phase  variable  measured  from  one  of  the  principal  axes  of  the  three- 
dimensional  trajectory.  If  and  only  if  the  negative  x-axis  aligns  with  the  projection  in  the 
x-y  plane  of  one  of  the  principal  axes,  we  can  equate  p  to  9  and  compare  Eqs.  (3.13), 
(3.14),  and  (3.17)  directly  to  Eqs.  (3.26). 

Likewise,  if  and  only  if  the  deputy’s  position  at  epoch  corresponds  to  one  of  the 
principal  axes,  we  can  equate  9  to  and  compare  Eqs.  (3.2)  directly  to  Eqs. 

(3.26) . 

However,  to  be  completely  general,  we  must  locate  the  principal  axes  of  the 
instantaneous  ellipse  in  three-dimensional  space  whose  existence  we  are  ultimately 
seeking  to  prove.  Again  following  Press’s  treatment  [71],  we  can  define  the  deputy’s 
distance  from  the  center  of  its  instantaneous  relative  trajectory  as 


Using  Eqs.  (3.13),  (3.14),  and  (3.17), 

d  (t)  = 
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If  the  instantaneous  relative  trajectory  is  indeed  an  ellipse,  then  d  (t)  will  have 
maximum  and  minimum  values  at  the  endpoints  of  this  ellipse’s  semi-major  and  semi¬ 
minor  axes,  respectively.  To  locate  the  extrema  of  d  (t) ,  we  can  set  its  p  -derivative 
equal  to  zero. 


d_ 

dp 


d(t)  =  0 


Differentiating,  and  assigning  P^  to  be  the  value  of  P  corresponding  to  the  1- 
direction  of  the  principal  axis  frame. 


ld{t) 


2a/sin/?,eos/?.-^sin/?,eos/?. 
+2z„ax'  sin  (/?,  +  r)  cos  (/?„  +  y) 


(3.27) 


Of  course,  because  Eq.  {3.21)  is  true  whenever  the  numerator  of  the  right-hand 
side  is  zero,  we  can  simplify  and  use  the  trigonometric  identities  in  Eqs.  (A. 2)  and  (A.  11) 
to  say 


0  =  3a/  sin cos 

-l4Zmax^  [sin  P^  cos  y  +  cos  P^  sin  y\ [cos  P^  cos  y  -  sin  P^  sin  y\ 


(3.28) 


Expanding,  rearranging,  and  applying  the  trigonometric  identities  from  Eqs.  (A.l) 
(A.  10),  and  (A.  1 1),  we  find  that 
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-Zmax"sin(2f) 


(3.29) 


tan(2;0j  =  - 

4«/+Zmax"cOs(2f) 

During  one  chief  orbit,  as  P  goes  from  zero  to  2n  ,  the  function  tan(2y9)  repeats 
four  times;  in  other  words,  tan(2y0)  is  periodic  with  a  period  of  ^  .  Thus  for  any  given 
value  of  tan  (2y0) ,  such  as  the  constant  right-hand  side  of  Eq.  (3.29),  the  solution  set  will 
contain  four  values  of  p  separated  by  ^  .  Since  we  have  defined  P^  to  correspond 
with  the  1 -direction  of  the  principal- axis  frame,  we  can  assign  the  second  solution 

to  correspond  with  the  2-direction.  The  remaining  two  solutions  will  correspond  with  the 
negative  1-  and  2- axes. 

Solving  Eq.  (3.29), 


-^max'sin 

^^«/+Zmax'cOs(2^)^ 

Press’s  analysis  [71]  stopped  after  finding  two  of  the  solutions  for  in  a  slightly 
different  form.  However,  we  must  continue  in  order  to  define  the  mapping  between 
reference  frames. 


=^tan 


50 


Since  our  goal  was  to  find  a  phase  variable  measured  from  the  principal- axis 
frame’s  1 -direction,  we  can  now  find  the  general  form: 

This  allows  us  to  re-write  the  deputy’s  c-frame  parametric  trajectory  equations 
(Eqs.  (3.13),  (3.14),  and  (3.17))  as  follows,  using  the  trigonometric  identities  from  Eqs. 
(A.2)and  (A.  11): 


X  ( t)  =  [cos  (3^  cos  6  -  sin  sin  -l- 

y  (t)  =  [sin  [3^  cos  0  -l-  cos  [3^  sin  6'\  +  (t) 

Z  (0  =  ^max  [sin  {/3^+r)cosO  +  cos{/3^+  y )  sin  O] 

This  form  can  be  compared  directly  to  Eqs.  (3.26),  showing  that  the  instantaneous 
relative  trajectory  viewed  in  three  dimensions  is,  in  fact,  an  ellipse,  with 


au  = 


-a„ 


2 


^maxSin(y0,+y) 


hv  = 


■  O 

ysin^^ 

a.  coS;5, 


2maxCOs(;5^+y) 


(3.30) 


(3.31) 
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yA^) 

0 

We  can  perform  two  simple  checks  on  these  results.  First,  we  could  find  au  and 
bv  by  substituting  the  solutions  for  and  (5^^  directly  into  Eqs.  (3.13),  (3.14),  and 
(3.17)  and  subtracting  the  coordinates  of  the  center  from  the  resulting  position 
coordinates;  the  results  match  Eqs.  (3.30)  and  (3.31). 

Second,  we  could  calculate  the  dot  product  of  au  and  hv  to  verify  orthogonality: 

au-bv  =  3a/  sin  cos  (3^ 

[sin  cos  Y  +  cos  sin  f  ]  [cos  P^  cos  y  -  sin  P^  sin  y\ 

Erom  Eq.  (3.28),  we  can  see  that  this  dot  product  is  indeed  zero. 

Note  that  we  now  have  a  general  method  for  transforming  from  principal- axis 
coordinates  into  c-frame  coordinates  and  vice  versa.  If  we  call  the  principal- axis  frame  r, 
the  necessary  rotation  matrices  are 

N  [dj 

R"'  =[«"' 

where 
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au 


u  = 


V  = 


\au\ 

bv 


Brief  Example. 

One  geometrical  insight  we  can  gain  using  the  results  of  this  chapter  concerns  a 
case  involving  only  in-plane  motion,  no  drift,  and  motion  centered  about  the  origin.  In 
this  case,  the  maximum  displacement  of  the  deputy  from  the  chief  equals  .  We  know 

(see  Eq.  (D.6))  that  the  deputy’s  orbit  radius  from  perigee  is  (l  -  ) .  Also,  since 

boundedness  indicates  that  the  chief  and  deputy  orbits  have  equal  semi-major  axes,  then 
we  can  say  the  distance  between  the  deputy’s  perigee  point  and  the  chiefs  orbit  radius 
equals  .  As  we  can  see  from  Figure  3,  this  perigee  displacement  equals  the  semi¬ 
minor  axis  of  the  relative  ellipse,  which  equals  ^ .  Finally,  we  have  an  expression  for 
the  maximum  distance  between  chief  and  deputy: 

Pmax  =«.  (3-32) 
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IV.  Parameterized  Models  of  the  Elliptical  Chief  Problem 


Objectives 

As  with  any  dynamics  model,  every  relative  satellite  motion  model  is  an 
approximation,  incorporating  certain  assumptions  that  affect  the  range  of  validity  for  the 
model.  For  example,  the  HCW  model's  assumption  of  a  circular  chief  orbit  means  that, 
for  a  given  application's  error  tolerance,  there  is  a  value  of  chief  orbit  eccentricity  beyond 
which  the  HCW  model  is  invalid. 

Each  model  represents  a  balance  between  accuracy  and  simplicity.  Choosing  the 
best  model  for  an  application  depends  on  getting  that  balance  right.  For  example,  if  a 
spacecraft  performing  RPO  is  to  have  an  autonomous  impulsive  maneuver  scheme,  its 
on-board  planning  computer  must  use  a  parameterized  relative  motion  model  that  can  be 
easily  optimized  (e.g.,  for  elapsed  time  or  velocity-change  magnitude)  to  determine  the 
timing,  magnitude,  and  direction  of  the  next  impulse.  Likewise,  if  the  proximity 
operations  mission  has  geometrically  defined  constraints,  such  as  field-of-view 
restrictions,  the  model's  parameters  should  have  a  simple  relationship  to  the  geometry  of 
the  relative  trajectory.  These  simplicity  requirements  must  be  balanced  against  accuracy; 
if  the  spacecraft  is  to  perform  RPO  near  satellites  whose  orbits  have  non-negligible 
eccentricity,  the  HCW  model  may  be  inadequate. 

Finding  a  good  relative  motion  model  whose  parameters  yield  significant 
geometric  insight  for  such  an  application  (an  impulsive  maneuver  scheme  for 
autonomous  RPO)  is  the  goal  of  the  present  chapter.  Having  seen  that  ROEs  constitute  a 
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satisfactory  model  for  the  circular  chief  problem,  we  will  examine  a  novel  method  (the 
Virtual  Chief  method)  to  see  whether  it  is  a  similarly  satisfactory  answer  to  the  elliptical 
chief  problem. 

In  other  words,  we  wish  to  relax  the  eccentricity  restriction  of  the  HCW  model. 
Specifically,  the  end  result  should  be  a  relative  motion  model  with  four  characteristics. 
First,  the  model  should  have  a  closed-form  solution  that  is  geometrically  intuitive, 
yielding  a  high  degree  of  insight  into  the  physical  trajectory  of  the  deputy  with  respect  to 
the  chief. 

Second,  the  model  should  lend  itself  to  a  high  degree  of  operational  efficacy.  In 
other  words,  developing  impulsive  maneuver  schemes  and  relative  navigation  algorithms 
should  not  involve  a  high  degree  of  complexity  or  computational  cost. 

Third,  assumptions  of  small  eccentricity,  small  separation  distances  (relative  to 
the  chief  orbit  radius),  and  Keplerian  two-body  dynamics  are  permissible.  Otherwise,  the 
model  should  not  exclude  any  classes  of  relative  orbit  problems.  For  example,  it  should 
accommodate  drifting  trajectories,  non-repeating  trajectories  which  experience  a  secular 
drift  due  to  a  difference  in  semi-major  axes.  (Of  course,  this  accommodation  may  last 
only  a  finite  time,  until  the  drift  causes  any  small-separation  assumption  to  be  exceeded.) 
Also,  the  model  should  permit  the  chief  and  deputy  orbits  to  have  widely  different 
arguments  of  perigee  and  right  ascensions  of  the  ascending  nodes  (if  the  orbits  are  near- 
equatorial),  provided  the  small-separation  assumption  is  not  violated. 

Fourth,  if  the  chiefs  eccentricity  is  zero,  the  model  should  reduce  to  the  HCW 
equations. 
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Possible  Approaches 


The  search  for  an  appropriate  relative  motion  model  could  have  proceeded  along 
any  of  five  possible  approaches. 

Neglect  Chief  Eccentricity. 

The  first  approach  would  be  to  use  the  HCW  model,  ignore  the  eccentricity  of  the 
chiefs  orbit,  and  accept  the  resulting  error.  This  approach,  although  used  for  some 
applications,  translates  into  fuel  penalties  even  for  small  eccentricities  [41]. 

Lawden  's  Equations. 

The  second  approach  would  be  to  use  Lawden's  equations,  attempting  to  find 
geometrically  meaningful  parameters  in  some  of  the  existing  solutions.  A  time-domain 
version  of  these  equations,  similar  to  the  notation  of  De  Vries  [26]  or  of  Inalhan, 
Tillerson,  and  How  [41],  can  be  written  as 
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Al  - 


^52 


/  \  3 

l  +  e^  cos 

\  J 


/  \  3 

1  +  cos  V(, 

V  ) 


(4.2) 


^3 


1  +  cos  V(, 

V  J 


As  a  linearized  two-body  model,  Lawden’s  equations  would  likely  yield  very 
accurate  results  over  small  relative  distances  and  short  timescales  up  to  moderately  large 
eccentricities  (i.e.,  around  0.6)  [58,  97].  However,  these  solutions  are  often  complicated, 
and  the  resulting  geometry  is  difficult  to  describe,  especially  if  drift  is  allowed. 

The  case  where  chief  and  deputy  orbits  have  the  same  semi-major  axis  (bounded 
motion)  is  a  different  matter.  Using  as  the  semilatus  rectum  of  the  chiefs  orbit, 

Sengupta  and  Vadali  [86]  gave  the  following  boundedness  condition  from  their  solution 
to  Lawden’s  equations  (converted  to  the  notation  of  this  thesis): 


0  =  — (2 
Pc 


-I-  cos  (t))  (l  -I-  cos  (0)"  ^(0  +  ^cJ—  sinvc 


— ^  srn  Vp  (t)  (l  + 
Pc 


^c  (Of  y(0  +  (^  +  ^c  cos  ^y(t) 

V  P 


After  evaluating  at  epoch  time  tg ,  applying  the  two-body  orbit  properties  in  Eqs.  (D.7) 
and  (D.l  1),  and  multiplying  by  constant  factor  ,  the  above  equation  becomes 


57 


2  +  erCO%VrQ  I  \2  Cr  (  2\X 

0  =  — ^ - r^(l  +  ecCOSVc,oj  ^o  +  ^sinVco(l-ec  )  K 


(i-v) 

■^sini'„  , 
(1-c^)  ^ 
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+  e^j  cos  Q  ^  Jo  “*" 


(4.3) 


-(l  +  ecCOSVc,o)  jo 


where  epoch  conditions  are  indicated  by  the  subscript  0. 

To  check  this  result,  we  can  use  an  alternate  but  equivalent  approach  from  the 
linearized  orbit-element  difference  method.  Bounded  motion  requires  simply  that  5a 
equal  zero.  Using  the  linearized  mapping  between  orbit-element  differences  and 
Cartesian  coordinates  from  Schaub  and  Junkins  [79],  after  some  manipulation. 


5a  = 


2(2  +  ecCOSVc(0) 
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Setting  5a  to  zero,  evaluating  at  ,  and  multiplying  by  the  constant  factor 


\ 

2 


(1 


causes  the  above  equation  to  equal  Eq.  (4.3),  showing  that  the  two  methods 


are  equivalent. 

As  we  have  seen,  for  bounded  motion  using  Law  den’s  equations,  there  has  been 
some  recent  success  describing  the  geometrical  impact  of  the  constants  of  integration  [54] 
and  defining  a  set  of  five  design  parameters  [86].  The  discovery  that  any  bounded- 
motion  solution  lies  on  a  definable  quadric  surface  [45]  is  certainly  geometrically 
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meaningful.  However,  because  an  RPO  maneuver  scheme  may  benefit  from  using 
drifting  trajectories,  these  options  were  not  pursued. 

First-  and  Second-order  Eccentricity  Approximations. 

The  third  approach  would  be  to  attempt  to  simplify  the  solution  to  Lawden's 
equations  by  making  a  further  approximation.  For  example,  a  series  expansion  of  the 
chiefs  true  anomaly  could  retain  only  terms  that  are  first-  or  second-order  in  the  chiefs 
eccentricity  (see  Eq.  (D.12)).  As  we  have  seen,  this  technique  has  been  used  for 
numerical  comparisons  often  in  the  literature,  and  Melton  applied  this  method  to  derive  a 
state  transition  matrix  of  relative  motion  which  can  be  initialized  at  any  time  during  the 
chiefs  orbit  [66] . 

Change  of  Reference  Frames. 

The  fourth  approach  would  be  to  attempt  to  achieve  a  simpler  representation  of 
the  relative  motion  by  using  some  chief-centered  reference  frame  other  than  the 
traditional  rotating  LVLH  frame.  For  some  applications,  a  variety  of  non-orthogonal 
reference  frames  have  proven  useful  [50].  Also,  using  an  inertial  coordinate  frame  (for 
example,  a  frame  parallel  to  the  chiefs  perifocal  frame)  has  shown  promise  in 
highlighting  simple  geometric  features  [56,  67],  but  this  idea  has  not  been  widely 
explored  in  the  literature. 

The  Virtual  Chief  Method. 

The  fifth  approach  is  the  Virtual  Chief  method.  It  is  an  attempt  to  avoid  the 
complexity  of  existing  elliptical-chief  models  and  take  advantage  of  the  simplicity  of  the 
HCW  model  and  its  closed-form  solution.  This  method  uses  a  virtual  chief,  a  third 
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satellite  (besides  the  chief  and  deputy)  whose  orbit  elements  are  the  same  as  those  of  the 
actual  chief,  except  for  eccentricity,  which  is  zero.  The  virtual  chief  is  a  theoretical 
object,  and  may  not  physically  exist.  It  is  defined  in  the  same  way  as  the  "reference 
platform"  of  Fasano  and  D'Errico's  higher-fidelity  model.  Their  model  is  first-order  in  the 
eccentricity  and  second-order  in  the  relative  distance,  and  it  incorporates  higher-order 
gravity  effects  [30].  The  linearized,  unperturbed  Virtual  Chief  method  is  the  approach 
selected  for  the  current  study,  although  adding  second-order  terms  could  constitute  an 
alternative  approach. 

The  Virtual  Chief  Equations 

The  Virtual  Chief  method  is  applicable  to  the  elliptical  chief  problem,  but  it 
assumes  a  small  eccentricity  for  the  orbits  of  both  chief  and  deputy.  This  eccentricity 
restriction  is  folded  into  the  small- separation  restriction,  because  we  assume  that  both  the 
chief  and  deputy  will  remain  near  the  virtual  chief.  Using  the  subscript  O  for  the  virtual 
chief  and  the  subscript  C  for  the  chief,  we  could  define  the  virtual  chiefs  orbit  in  terms  of 
semi-major  axis  a  ,  inclination  i ,  eccentricity  e ,  right  ascension  of  the  ascending  node 
Q ,  and  mean  anomaly  M  ,  as 
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Note  that  even  though  there  is  no  perigee  for  the  circular  orbit,  we  have  defined 
the  mean  anomaly  to  be  equal  to  the  chiefs  mean  anomaly.  Thus,  at  the  chiefs 

perigee,  both  mean  anomalies  will  equal  zero. 

Reference  Frames  for  the  Actual  and  Virtual  Chief 

Let  there  be  an  LVLH  reference  frame  c  fixed  to  the  chief,  and  an  LVLH  frame  o 
at  the  virtual  chief.  By  the  definition  of  the  virtual  chief,  its  orbit  is  coplanar  with  that  of 
the  chief.  Thus, 


=  o 


3 


As  a  result,  the  angular  separation  between  c  and  o  can  be  described  by  a  rotation 
about  the  third  axis  through  an  angle  Am  ,  representing  the  difference  in  arguments  of 
latitude  between  the  chief  and  the  virtual  chief  (see  Figure  4). 

Since  both  arguments  of  latitude  are  measured  from  the  ascending  node  shared  by 
both  orbits,  and  since  and  (the  chiefs  true  anomaly)  are  both  measured  from  the 
chiefs  perigee  point,  we  can  say  that 


Uq  —  COf.  +  Mg 

Ug=COg+  Vg  (4.5) 

Au  =Ug  -Ug  =  Vg  -Mg 
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Finally,  since  we  know  that  both  mean  anomalies  are  always  equal,  we  can 
construct  the  following  rotation  matrix  that  will  convert  o-frame  coordinates  into  c-frame 
coordinates: 


RO^C 


cos(v^-M(.)  sin(Vp-M^)  0 

-sin(Vp-M(,)  cos(v^  -M^)  0 
0  0  1 


62 


Let  the  angular  velocity  of  the  c-frame  with  respect  to  the  o-frame  be  co, .  Since 

/o 

both  orbits  will  always  be  coplanar,  we  know  that  the  only  component  of  a,  will  be  in 

A? 

the  out-of-plane  direction.  In  fact,  the  magnitude  of  this  component  will  be  the  time  rate 
of  change  of  Am  : 


CO,  = 


fd  ^ 
—  Am 

<^3  = 

S 

<] 

\dt  ) 

ydt  ) 

Finding  the  chiefs  mean  anomaly  from  Eq.  (D.3),  and  using  as  the  time  derivative  of 
the  chiefs  true  anomaly,  we  can  calculate  the  time  rate  of  change  of  Am  as  follows: 

dt  dt^  ^  ""  dA  pJ)  c  c 


Thus, 


The  angular  velocity  of  the  c-frame  with  respect  to  inertial  space  is  simply 


Taking  an  inertial  derivative  and  using  A  the  second  time  derivative  of  the  chiefs  true 
anomaly  yields  the  angular  acceleration: 
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(4.6) 


d  ^ 


If  Hq  is  the  virtual  chiefs  mean  motion,  then  examining  the  o-frame  in  a  similar  way 
results  in 


Likewise, 


d  ^  n 


(4.7) 


Vector  Equations. 

In  order  to  write  the  vector  equation  of  relative  motion  for  our  model,  we  must 

find  an  expression  for  — ^  r„/  .  If  the  positions  of  the  chief  and  deputy,  respectively, 

dt  /c 

from  the  virtual  chief  are  f^,  and  ,  then  we  can  also  make  the  simple  vector 

Vo  /o 


subtraction  statement  that 


Taking  an  o-frame  time  derivative, 


(4.8) 
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d  ^  "  d  ^ 

j/%  J/% 


d  ^ 


A  further  o-frame  time  derivative  yields 


dt^  %  dt^  % 


X 


Rearranging  will  provide  the  nonlinear  vector  equation  of  relative  motion  we  have 


been  seeking: 


c 


d^  ^  "  d  ^ 

— -  r^/  -  2co,  X  —  r„/ 

dt^  %  /o  dt  /c 


xi^/  (4.9) 

/c 


Eq.  (4.9)  is  the  Virtual  Chief  vector  equation  of  motion.  It  now  remains  to 
provide  scalar  values  for  the  c-frame  components  of  each  term. 

Scalar  Equations. 

First,  let  us  denote  the  position  coordinates  of  the  chief  and  deputy,  respectively, 
with  the  capital  subscripts  C  and  D;  and  let  us  do  likewise  for  the  o-frame  relative 

°  d^ 

velocity  and  acceleration  components  (for  example,  — ^  will  be  written  as 

dt  /o 

Next,  let  us  make  use  of  the  HCW  equations  (Eqs.  (3.1))  to  describe  the  motion  of 
the  chief  and  deputy  relative  to  the  virtual  chief.  Since  the  virtual  chief  is  in  a  circular 
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orbit,  the  assumptions  of  the  HCW  equations  are  not  violated,  as  long  as  the  magnitudes 


of  fr/  and  FU/  remain  small.  Thus,  we  find  that 
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Of  course,  since  the  orbits  of  the  chief  and  the  virtual  chief  have  the  same  semi¬ 
major  axis,  it  follows  that 


Furthermore,  since  we  know  that  the  orbits  of  the  chief  and  the  virtual  chief  are  coplanar, 
we  can  say  that 


Zc=Zc=Zc=0 


In  addition,  since  we  know  — ^  Fj,/  and  — ^  only  in  o-frame  coordinates,  we  must 

dt  /o  dr  /o 


apply  the  rotation  matrix  before  substituting  into  Eq.  (4.9).  In  other  words. 
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Next  we  must  find 


—  a, .  We  can  again  make  a  simple  vector  subtraction  statement, 
dt  /o 


Taking  a  c-frame  time  derivative, 


Substituting  from  Eqs.  (4.6)  and  (4.7), 


d  ^ 
dt  /o 


0 

0 


Now  we  are  prepared  to  rewrite  Eq.  (4.9)  in  matrix  format.  Substituting  from 
results  above. 
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Collecting  terms  and  performing  the  cross  products. 
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i* 
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We  can  derive  expressions  in  terms  of  x,  y,  and  z  to  replace  the  factors  ) , 

(xq  -  ip ) ,  ( jo  -  jc )  >  and  in  Eq.  (4. 10)  by  using  the  rotation  matrix  .  From  Eq. 
(4.8), 
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Since  R"^‘^  is  a  rotation  matrix,  and  therefore  orthogonal. 
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Inserting  the  appropriate  component  values  and  evaluating. 
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Next,  we  rotate  the  o-frame  derivative  of  Eq.  (4.8)  in  the  same  way: 
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(4.13) 
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Inserting  the  appropriate  values, 
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Calculating, 


X^-Xc 

cos(vc -Mc)[x-{vc-nc)y']-^m{vc -Mc)[y  +  {vc-n^)x'\ 
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_^D~  _ 

o 

z 

(4.14) 


Now  we  can  substitute  from  Eqs.  (4.12)  and  (4.14)  into  Eq.  (4.10): 


2(^c- 

-nc)y  +  {vc-ncf  x  +  v^y 

cos(V(, 

sin(V(^  -M^) 

O' 

y 

= 

-2(vc 

-nc)x  +  {vc-nc)^  y-Vc^ 

+ 

-sin(V(^  -Mf,) 

cos(Vc  -M^) 

0 

J'J 

c 

0 

0 

0 

1 

\c 


2nc{ 

sin(vc-Mc)[ 

x-{vc-nc)y~] 

+3«c^ { 

+  cos(vc-3^c) 

cos(V(,  -Mc)x- 

[y+i 

-sin( 

S' c  ^c) 

c 

4) 

)4 

-In^ 

{cos(vc-3^c! 

)\_x-{vc-nc)y 

']-%m{vc-Mc) 

l[j  + 

Sc->^c] 

Id) 

-tic 

2 

Z 
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Simplifying  and  performing  the  matrix  product  gives 


X 


y 


z 


\c 


2Vf.y  +  (^Vf.-n^y  x  +  V(,y  +  2n^  -n^'jx  +  cos^ )x 

-3n^^  cos(v^  -M^)sin(vp  -M^^y 

-2VcX  +  {vc-ncf  y-  v^x  +  2nc  (v^  -  )  3^ 

-3?i^^sin(Vp  -M^)cos(v^  -M^'^x  +  3n^^  sin^(v^  -M^)y 

. 2 . 

-He  Z 


(4.15) 


This  is  a  set  of  scalar  equations  of  relative  motion  that  depend  on  time  (through 
Me,  Vp  ,  ,  and  )  and  known  constants  and  .  Note  that  the  out-of-plane 

equation  is  identical  to  the  HCW  out-of-plane  equation.  Eq.  (4.15)  can  be  further 
simplified  and  written  conveniently  as  a  linear  time- varying  system  in  the  following  form: 


■  0 

0 

0 

1 

0 

o' 

X 

y 

0 

0 

0 

0 

1 

0 

y 

z 

0 

0 

0 

0 

0 

1 

z 

X 

A41 

^42 

0 

0 

2vp 

0 

X 

y 

Al 

^2 

0 

-2vp 

0 

0 

y 

_2'j 

c 

0 

0 

2 

-He 

0 

0 

0 

z 

L  Jc 

(4.16) 


where 


A4j  =Vc^  ~nc  +3n(?  cos^iv^ 

A42  =v^  -3n^^cos(Vp  -M ^)sm{y^  ~^c) 
A^j  =  -i/p  -3nJ'  sin(Vp  -M^)cos{v^  ~^c) 
A52  =  Vq  -nj'  +3nJ'  sin^(Vp  -M^^) 
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The  time  dependence  of  the  system  in  Eq.  (4.16)  is  implicit.  The  time  dependence 
through  could  easily  be  made  explicit  through  Eq.  (D.3).  However,  showing  the 

time  dependence  through  and  its  derivatives  requires  solving  Kepler’s  transcendental 
equation  (Eq.  (D.4)). 

Note  that  if  we  wished  to  account  for  differential  drag,  thrust,  or  other  forces 
affecting  the  relative  motion,  we  could  do  so  by  adding  the  appropriate  terms  to  the  right- 
hand  side  of  Eq.  (4. 16).  Note  also  that  this  model  resembles  the  time-domain  Lawden’s 
equations,  Eqs.  (4.1)  and  (4.2),  but  with  additional  terms. 

Einally,  note  that  Eq.  (4.16)  reduces  to  the  HCW  equations  when  the  chiefs 
eccentricity  is  zero.  To  show  this,  simply  substitute  the  following  circular-chief 
equivalent  values  for  the  chief  true  anomaly  and  its  successive  time  derivatives:  , 

and  0. 

The  Virtual  Chief  Solution 

Eq.  (4.16)  represents  a  system  of  second-order  differential  equations  whose 
dynamics  matrix  A(t)  is  known.  When  the  closed-form  solution  to  the  system  takes  the 

form  of  a  state  transition  matrix  d)  (t,  tg ) ,  that  matrix  must  obey  the  following: 

6(t,tg)=A(t)0(t,tg)  (4.17) 

The  initial  condition  for  Eq.  (4.17)  is  that  is  the  identity  matrix  [73]. 
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However,  rather  than  solve  the  differential  equations  directly,  we  shall  once  again 
take  advantage  of  the  known  HCW  motion.  The  HCW  position  and  velocity  solutions, 
whose  state  transition  matrix  forms  are  shown  in  Eqs.  (3.3)  and  (3.5),  can  be  applied  to 
the  o-frame  motion  of  the  chief  relative  to  the  virtual  chief: 


A  similar  construction  can  be  used  to  define  the  motion  of  the  deputy  relative  to  the 
virtual  chief. 


Vector  Solution. 

To  find  the  solution  in  vector  form,  we  begin  by  differencing  the  chief  and  deputy 
motions  expressed  in  Eqs.  (4.18)  and  (4.19)  as  follows: 


(4.20) 
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"  d  ^ 
dt  /o 


^  d  r 

(0- 


(4.21) 


For  initial  conditions,  we  can  evaluate  Eqs.  (4.11)  and  (4.13)  at  epoch  time  : 


{h)-rc/{h)  =[^'’"'(^o)7  r^/{h) 

'  /  O  Jo  L  /C  Jc 


dt  7c 


(4.22) 


Note  that  each  side  of  the  first  of  Eqs.  (4.22)  equals  rj^/  [tJ  ,  and  each  side  of  the 

L  /c  Jo 


second  of  Eqs.  (4.22)  equals  — (tg)  .  Substituting  into  Eqs.  (4.20)  and  (4.21), 

dt  /c 


o;  'D/ V'O 


T  ^  d 

■*'®rv  (^o)]  (^o)^  0^  (^o) 


(4.23) 


7%7-  7%7\ 


D/^VO 


(  ) 


(4.24) 


+®c/(^o)xv(^o) 
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Then  we  can  substitute  Eqs.  (4.23)  and  (4.24)  into  rearranged  forms  of  Eqs.  (4.11) 


and  (4.13),  yielding 

%(^o)  ^ 

Y —  r 

+  (0]  ®  ^  ^0  )  (^0  )]^ 

^+d5^(to)xy^(to) 

-ic 

T— r  (n 

+  [r‘’-^  (0]o„  (t„)]"  dt  Vc  (4.26) 

+®c/(^o)xV  (^o) 

/o  /C  Jc 

-  ®./(0^v(0 

_  /o  /C  Jc 


Einally,  substituting  Eq.  (4.25)  into  Eq.  (4.26), 
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Eqs.  (4.25)  and  (4.27)  are  the  solution  we  have  been  seeking,  in  vector  form.  The 
relative  position  and  velocity  at  any  time  are  given  as  a  function  of  initial  position,  initial 
velocity,  time,  and  a  variety  of  time-varying  elements  from  the  rotation  matrices  and 
angular  velocity  vectors. 

Scalar  Solution. 

The  elements  of  all  the  matrices  and  vectors  in  Eqs.  (4.25)  and  (4.27)  are  already 
known.  Eluding  the  scalar  solution  equations  is  then  a  straightforward  matter  of  vector- 
matrix  operations.  Eq.  (4.25)  becomes 


x{t) 

cos[v^{t)-M^{t)]  sin[v,( 

0-^c(0] 

o' 

yit) 

-sin[v,(?)-M^(0]  cos[v,( 

0-^c(0] 

0 

0 

0 

1 

-3 cos[n 

c(^- 

-0]  +  4 

0 

0 

6sin[«^  (r- 

0] 

)  1 

0 

0 

0 

cos  [t  - 

cos[v,„-M 

c,o  ] 

-sinKo 

-^c,] 

o' 

^0 

sin[v,„-M 

CO  ] 

cos[v,„- 

0 

>'o 

0 

0 

1 

cos[vc{t)-M^{t)]  sin[v,(0-M^(0]  0 

-sin[v,(r)-M^(?)]  cos[v^{t)-M^{t)]  0 

0  0  1 

sin (?-?„)]  2cos[n^{t-tJ]  2 


2 cos  (?-?„)]  2  4  sin  (?-?„)] 


-3{t-0  0 


COS [v^ , 

-sinfv^„  1 

o' 

[41 

0 

L  c,o 

C.O  J 

L  C,0  C,0  J 

0 

0 

sin[v,„-M^J 

cos[v,„-M^J 

0 

4 

+ 

0 

X 

> 

0 

0 

1 

1 

_4,o-«c. 

_  '^0  _ 

(4.28) 
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For  the  sake  of  brevity  in  the  following  development,  we  introduce  constant 
intermediate  parameters  A,  B  ,  C  ,  and  D  such  that  the  initial  conditions  can  be 
expressed  as 


A 

B 


1 

1 

O 

1 _ 

= 

ticD 

ticC 

_1 

o 

_  ^0  _ 

We  can  find  the  values  of  the  intermediate  parameters  by  calculating  Eq.  (4.22): 


C  =  —i 
n. 


A  =  cos  [vc.o  -  ^c,o  ]  ^0  -  sin  [ -  M ^  ]  y„ 

B  =  sin  [vc,o  -  o  ]  ^0  +  cos  [ v^.o  “  ]  3^o 

[ sin [vc,o  - ] [io  )  Jo ] 

[^C,0  ~^C,o\  jo  ■*“('^C,0  ~^c)^o 
cosj^Vp Q  —  Q  J  Vo  —  (vp.o  Jo] 
[^C,0  ~^C,0~\  Jo  "^(^C.o  ~^c)^Q 


D  =  —\ 
n. 


-sm 


(4.29) 


Using  these  parameters  to  evaluate  Eq.  (4.28),  we  find  that  the  three  position 
components  are 
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^(0  =  sin[vc(0-^c(0]l 


+  COS 


\vc{t)-Mc{t)\\ 


|6sin[np  -6n(,  A  +  fi 

+  {4sin[nc(?-?o)]-3nc(?-?o)}c 
+  |2cos[nc 

|-3cos[n(,  (?-?o)]  +  4|  A 
+  |-2cos  [n^  {t  -  tfj )]  +  2|  C  +  sin  [n^  [t  -  )]  D 


(4.30) 


y(?)  = -sin[V(.  (?)  •< 


+  COS 


[vc{t)-M^{t)\\ 


|-3cos  [n^  {t  -  ?o )]  +  4|  A 
+ 1-2  cos  [n^  (?  -  ?o )]  +  2|  C  +  sin  [n^  (?  -  ?o )]  -D 
|6  sin  (?-?„)]-  (?-?g)|A  +  5j 

+  {4sin[nc(?-?o)]-3nc(?-?o)}c 
+  |2cos[nc  (?-?o)]-2|d 


(4.31) 


/  N  r  /  N-|  Sin (?-?„)] 

z(?)  =  cos[nc(?-?o)]zo+ - ^ - ^4 


(4.32) 


The  velocity  components  are  a  little  more  complicated.  Eq.  (4.27)  becomes 
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■.(4] 

1 

o 

o 

n 

1 

sin[v^  (»)]  0“ 

y(0 

= 

-sin[v^(r)-M^(f)] 

cos[v^  (»)]  0 

_i(oJ 

0 

0  1_ 

3”c  ®“[”c  ('-'o)] 

0 

0 

COS  V  -  M 

>-  C.O  C,0  J 

-sinFv  -M  1  0 

>-  C.O  C,0  J 

0 

6«^  cos[n^  ('  ~  )]  ^ 

0 

0 

sin  V  -  M 

>-  C.O  C,0  J 

COS  V  -  M  0 

>-  C.O  C,0  J 

0 

0 

-«c 

0 

0  1_ 

-^0- 

cos  (0-M^(0]  sin  (0-M^(0]  0l 
+  -  sin  [v^  (0-M^(0]  cos  (0-M^(0]  0 

0  0  1_ 
cos[w^  (f-r^)]  2sin  [rt^  (r  - /^ )]  0 

— 2  sin  [  ( r  -  ?^  }  ]  4  cos  —  0 

0  0  cos  [rt  - 1  )] 


cos  V  -M  -sin  v  -  M  0 

>-  C,0  C.O  J  >-  C,0  C,0  J 


sin  V  -  M 

>-  C,0  C.O  J 

0 


4v 


-M  1  0 

,0  C.O  J 

0  1 


cos[r^  (f)]  sin[r^  (r)-M^  (f)]  0 

-sin[r^  (r)-M^  (f)]  cos[r^  (r) (f)]  0 

0  0  1 

-3cos[n  (r-?)]  +  4  0  0 

6sinr«0  — 1  0 

Lc''  0-'-'  C^  0-' 

0  0  cos  [n  (f  —  f  }] 


V  fO 


cos  \v  -  M 

L  C.O  C.O  J 

-sinFv  -M  1  0 

L  C.O  C,0  J 

X 

0 

sin  V  -  M 

L  C.O  C,0  J 

COS  \  V  -  M  0 

L  C.O  C,0  J 

0 

0  1_ 

-=^o- 

cos  (0-M^(0]  s  n  [v^  (0-M^(0] 

+  -sin[F^  (/)]  cos[i/^  (0-M^  (0] 

0  0 


0 

0 

1 


v'vy 


sin  Fn^  ( f  —  f 

„)]  -2cos[n^(f 

2  COS  Frt  C  ?  -  ^ 

)]-2  4sin[n^(r-f^ 

0 

0 

cos  \  V  -  M 

-  sin  V  -  M 

L  C.O  C.O 

J  L  C.O  C,0  J 

sin  V  -  M 

COS  V  -  M 

>-  C.O  C,0 

-•  >-  C.O  C,0  J 

0 

0 

0 

0 


sin  [n^  (r 

- 1 

0 

)]J 

0 

0 

0 

X 

V  -  n 

•—  C,0  C  — 1 

-^0- 
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Evaluating  with  A,  B  ,  C  ,  and  D , 


x(t)-sin 


btif,  cos  [n^  (t  -  to )]  -  6n^  |  A 

|4cos  [«(,  (t  -  to )]  -  3|  C  -  2no’  sin  [n^.  (t  -  to )]  D  J 


+  COS 


-sin 


Slip  sin[no-(t-to)]}^ 

\^+2n^  sin  (t  -  to )]  C  +  cos  [n^;  (t  -  to  )]o 
|-3cos  [lie  (^  “  h )]  +  A 
+  |-2cos[no-(t-to)]  +  2|c 


in[vc(t)-Mc(t)_ 


+  COS 


Vc{t)-Mc{t)_ 


+  sin[n^  (t-to)]D 

|6sin[n^  (t-to)]-6«^  (t-to)|  A  +  5 

+  {4sin  [n^  (t  -  to )]  -  3no;  (t  -  to )}  C 
+  |2cos  [wo-  (t  -  to )]  “  ^ 


=  -  sin 


+  COS 


t'clt)  ■^c(t)_ 


-{’^c(t)-«c) 


sin 


{3«cSin[«c(t-to)]}A 
+2«c  sin  \ji(,  (t  -  to )]  C  +  lie  ‘tos  [iic  (t  “  h 
jdile  cos  [lie  (t“to)]-6ile}  ^ 

+iie  |4cos  [lie  (t  “  to )]  -  3|  C  -  2n^  sin  [iie  (t  “  tg )] 
|6  sin  [lie  (t  “  to )]  “  6iie  (t  “  tg  )|  A  +  5 
+  {4  sin  [lie  (t  -  to )]  -  3iie  (t  -  to )}  C 
{2  cos  [lie  (t- to)]- 2}  D 
|-3cos[iie(t-to)]  +  4|  A 


D 


_t'c(t)  ■^c(t)_ 


+  COS 


_t'c(t)  ■^c(t)_ 


+  |-2cos[iie(t-to)]  +  2|c 

[iie(t-to)]Z) 


+  sin 


(4.33) 


(4.34) 
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z  (? )  =  sin  (?  -  ?o  )]  ^0  +  cos \nc{t- )]  ^ 


(4.35) 


Thus,  Eqs.  (4.30)-(4.35)  are  the  Virtual  Chief  solution.  A  comparison  reveals 
that  Eqs.  (4.33)-(4.35)  are  indeed  the  time  derivatives  of  Eqs.  (4.30)-(4.32). 

Note  that  the  out-of-plane  solution  is  identical  to  HCW,  and  that  if  the  chiefs 
eccentricity  is  zero,  Eqs.  (4.30)  and  (4.31)  reduce  to  the  in-plane  components  of  Eqs. 
(3.2). 

State  Transition  Matrix. 

We  can  arrange  the  Virtual  Chief  solution  into  a  new  state  transition  matrix 
d)(t,to) ,  such  that 


x{t) 

Xg 

y{t) 

3^0 

z{t) 

II 

3^ 

Zo 

x{t) 

Xg 

y^) 

^0 

_z{t) 

c 

Jo_ 

The  matrix  is 


O(t,to) 


®12  0  ®14 

®21  ®22  0  O24 

0  0  cos[nc  (t-tg)]  0 

®41  ®42  ^  ®44 

®51  ®52  0  O,, 

0  0  -fig.  sin[?i(^  (t-tg)]  0 


®15  0 

®25  0 

Sin[flg-(t-tg)] 


®45  0 
0 

0  COs[fIg,  (t-tg)] 


(4.36) 
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The  sixteen  in-plane  elements  of  O  (tTo)  would  be  too  cumbersome  to  list  here; 

they  are  found  in  Appendix  E.  It  can  be  shown  that  d)  (t,  tg )  satisfies  Eq.  (4. 17),  and  that 

the  time  derivatives  of  Eqs.  (4.30)-(4.35)  are  consistent  with  the  differential  equations, 

Eq.  (4.16).  Thus,  the  solution  presented  here  is  the  unique  solution  to  the  linear  time- 
varying  system  of  the  previous  section  [24] . 

Parameterization. 

Our  goal  is  to  find  a  meaningfully  parameterized  model.  Since  the  Virtual  Chief 
out-of-plane  motion  is  identical  to  HCW  out-of-plane  motion,  we  can  adopt  two  ROEs, 
Zmax  ¥o  (defined  in  Eqs.  (3.15)  and  (3.16)),  and  re-write  Eq.  (4.32)  as 

z(0  =  ^maxSin(nc(^-^o)  +  n)  (4-37) 


We  can  rewrite  the  in-plane  position  equations  in  the  following  manner,  using  the 
constants  of  Eq.  (4.29),  as  well  as  the  trigonometric  identities  in  Eqs.  (A. 7)  through 
(A.  10): 


KO  =  (^|^jsin[^c(0  +  Vp-Vo]  +  (^-^^-3cjcos[vc(0  +  V;,-Vo] 

-  2D)  sin  (t)  -  n^t  +  ]  -l-  (4A  -l-  2C)  cos  (t)  -  n^t  +  ] 


v2  j 


sin  +  - 


-A  +  C 
2 


cos 


/c(0“2v  +  Vp+Vo_ 


-(6A-l-3C)np  (t-tg)sin  +  n^t ^ 
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j(0  = 


(  9 
-A  +  3C 

U 


3 


sin[vc(0  +  Vp-Vo]  + 


flol 
U  J 


cos 


7c(0  +  Vp-Vo_ 


-(4A  +  2C)sin  V(.{t)-n^t  +  n(.tp  +(5-2D)cos  v^{t)-n^t  +  n^tp 


^  3  \ 

-  A  +  C  J  sin  [vp  (? )  -  2nc?  +  ]  +  -  D  cos  [v^  (? )  -  2n, 

(6A  +  3C)n^  (?-?q)cos[vp  +  J 


-2V  +  Vp+Vo_ 


Next,  we  can  apply  the  Harmonic  Addition  Theorem  (see  Appendix  B),  yielding 


: (t)  -  D V  ^  A V  27 AC  +  9C"  sin 


-  4BD  +  4D"  +  16A"  +  16AC  +  4C"  sin 


(1  9 

'4  4' 


+J-D"+-A"  +  3AC  +  C"  sin 


^c(0  +  Vp-Vo 
+aton2(-3A-  2C,D) 

Vc(t)-V  +  Vp 
+aton2  (4A  +  2C,  B  -  2D) 

Vc(t)-2v  +  V/,  + Vc 
+aton2(3A  +  2C,D) 


(4.38) 


-(6A  +  3C)np  (t-to)sin  + 


/-  D  V  —  A  V  27 AC  +  9C"  cos 

'^c(0  +  V/,-Vo 

H  4 

+aton2(-3A  -  2C,  D) 

+^JB^  -  4BD  +  4D"  +  16A"  +  16AC  +  4C"  cos 


1  9 

+,  -DV-A"  +  3AC  +  C"  cos 
\4  4 


Vc{t)-nct  +  nctp 

+aton2(4A  +  2C,  B  -  2D) 

^c(0“2V  +  Vp+Vd 
+aton2(3A  +  2C,D) 


(4.39) 


-(6A  +  3C)n(^  (t-tQ)cos  n^t  + 


Let  us  define  the  following  parameters: 
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(4.40) 


A,=yD^  +  {iA  +  2C)'‘ 

A  =  +  atan2(^3A  +  2C ,  D) 

A^  =  ^{B-  2Df  +  (4A  +  2cf 
^2  =  atan2(^4A  +  2C,B-  2D) 
A3  =  -(6A  +  3C)np 


Obviously,  these  five  parameters  are  not  independent,  but  are  defined  by  the  four 
intermediate  constants  A ,  B  ,  C  ,  and  D ,  which  are  defined  by  the  four  in-plane 
Cartesian  initial  conditions.  Substituting  these  parameters  into  Eqs.  (4.38)  and  (4.39), 
and  using  the  definition  of  mean  anomaly  in  Eq.  (D.3),  as  well  as  the  identities  in  Eqs. 
(A.  12)  through  (A.  15), 


x(t)-AiSin[vc(0  +  Vp- 

+  ^A^sm[vc{t)-2nc 


y(t)-AiCOs[vc(0  +  Vp- 

+  \\^o%\vc{t)-2nc 


>1  ]  +  4  sin  [  Vc  (0  “  +  ^ctp  +<l>2~] 

^  +  V/,+^^i]  +  4(^-^o)sin[vc(0“V  +  Vj 

-  ]  +  A2  cos  [vc  {t)  -  n^t  +  +(l)^~\ 

^  +  Vp+^^i]  +  4(^-^o)cos[vc(t)-V  +  Vj 


(4.41) 


(4.42) 


The  parameterized  solution  as  it  now  stands  comprises  Eqs.  (4.41),  (4.42),  and 
(4.37).  It  requires  knowledge  of  chief  mean  motion  and  time  of  chief  perigee  passage. 
(Of  course,  calculating  the  parameters  from  Cartesian  initial  conditions  also  requires 
knowledge  of  chief  true  anomaly  at  epoch  and  chief  angular  rate  at  epoch.)  The  seven 
parameters  ( A, ,  (1)^,  A^,  ^2  >  4  ’  ^max  ’  reduced  to  six  by  defining  A3  in 
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terms  of  the  others.  The  simplest  method  for  doing  so  is  to  examine  the  right  triangle 
suggested  by  the  definition  of  .  Its  hypotenuse  would  equal  ,  and  sin  (j)^  would  equal 

(4A  +  2C)/4.  Thus, 


—  A2  sin  ^2  =  6A  +  3C 


(4.43) 


We  ean  substitute  Eq.  (4.43)  direetly  into  the  definition  of  ,  yielding 


We  ean  reeognize  some  of  the  terms  in  Eqs.  (4.41)  and  (4.42)  as  representing  the 
differenee  in  argument  of  latitude  between  the  ehief  and  the  virtual  ehief,  Am  (t) .  Erom 
Eqs.  (4.4),  (4.5),  and  (D.3), 

Vp  (t)  -  n^t  +  =  Am  (t)  (4.44) 

This  allows  us  to  define  the  relative  motion  with  six  Virtual  Chief  parameters.  A, , 
(j)^,  Aj ,  ^2^2’  ^max  ’  ¥0  ■  The  position  equations  beeome 
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(4.45) 


x{t)  =  Aj  sin  [ Am  (?)  +  n^t-(l)^~^  +  sin  [Am  (?)  +  ^2^2] 

1  3 

+  -  Aj  sin  [Am  (?)  -  Mp?  +  ^2^1  ]  -  —  A2  sin  (l)^n(.  (?  -  ?o )  sin  [Am  (?)] 

y  (?)  =  Aj  cos  [  Am  (?)  +  Mp?  -  ]  +  A2  cos  [Am  (?)  +  ^2  ] 

1  3 

+  -  Aj  cos  [ Am  (?)  -  Mp?  +  ^2^1  ]  -  —  A2  sin  (l)^n^  (^  “  (0] 

z(0  =  ^maxSin(nc(?-?o)  +  n) 


Thus,  the  in-plane  motion  is  seen  to  operate  in  three  frequencies  and  involve 
mixed-secular  terms  of  the  form  sin  («2^  -t-  «3 ) .  This  form  is  similar  to  the  terms  that 

concerned  Melton  in  his  own  time-explicit  solution  [66].  However,  in  Eqs.  (4.45)  it  is 
nothing  more  than  the  projection  of  the  deputy’s  drift  (which  is  along-track  only  in  the 
virtual  chiefs  frame)  into  the  actual  chiefs  frame. 

The  Virtual  Chief  parameters  can  be  used  to  write  the  relative  velocity  solution,  as 
well.  Taking  a  time  derivative  of  Eqs.  (4.45),  using  the  Harmonic  Addition  Theorem  and 
the  trigonometric  identities  from  Eqs.  (A. 7)  through  (A.  10),  and  then  rearranging  terms 
yields 


v(?)  =  -A2M^yl-  — sin^  (j)^  cos[AM(?)-f  a?aM2(-sin^^2’2cos(2i2)] 

2  3 

-  -  AjMc  cos  [  Am  (?)  -  Me?  +  ]  +  -  A2  sin  (?  -  ?o )  cos  [  Am  (?)] 

-VVf,  (?)  Aj  cos  [Am  (?)  -t-  Me?  -  ]  -I-  fe  (?)  A2  cos  [Am  (?)  +  ^^2  ] 

1  3 

+  3^c(OAcos[AM(?)-Me?  +  ^2il]--fe(04sin^2«e(^-^o)cos[AM(?)] 
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(4.47) 


y(t)  =  A2n(,  .^1  -  ^  sin  ^  ^2^2  sin  [  Am  ( ? )  +  a  ton  2  ( -  sin  ^2^2 , 2  cos  ^2^2)] 

2  3 

+  —  sin  [ Am  (?)  -  m^?  +  sin  (?  -  ?o )  sin  [ Am  (? )] 

-Vp  (0  A  sin[AM(?)  +  M(,?-^zii]- A  (0  A  sin[AM(?)  +  ^zi2] 

1  3 

-3'^c(0Asin[AM(?)-Mc?  +  ^^l]  +  -A(0AsinA«c(^-^o)sin[AM(?)] 


2  (0  =  ^max'^c  cos  (m^  (?  -  ?0  )  +  ) 


(4.48) 


Eqs.  (4.46)-(4.48)  can  also  be  derived  direetly  from  Eqs.  (4.33)-(4.35). 

The  set  of  Virtual  Chief  parameters  is  a  realization  equivalent  to  the  six  Cartesian 
initial  conditions,  and  there  is  a  one-to-one  mapping  between  the  two  sets.  Using  the 
ehief  s  orbit  elements  as  given  values,  the  Virtual  Chief  parameters  ean  be  ealculated 
from  Cartesian  initial  conditions  via  Eqs.  (3.15),  (3.16),  (4.29),  and  (4.40).  The  Cartesian 
initial  conditions  can  be  caleulated  from  the  Virtual  Chief  parameters  by  evaluating  Eqs. 
(4.44).(4.48)  at  t, : 


AMq  V^  g  M(2?q  -I- 


Xg  =  A,  sin  [  Amq  -I-  n^tg  -  ]  -l-  A2  sin  [Am,,  -l-  ^2^2  ]  +  ^  A^  sin  [  Amq  -  n^-tg  +  (j)^  ] 


yp  =  Aj  cos  [  Amq  -I-  n^tg  -  ]  -l-  A2  cos  [  Amq  -l-  ^2  ]  +  ^  A  cos  [  Amq  -  n^tg  +  ] 


z^^  =  z^^my/g 
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/  3  2 

X,)  =  -A^n^Jl-^sin^  (j)^  cos[AMg  +aton2(-sin^2’2cos^zi2)]“  — A^c  cos[Amo  -n^t^  +  ^zij] 

+^c,o  A  cos  [  Amq  +  n^tQ  -  ]  +  A,o  A  cos  [  Am^  +  A  ]  +  ^  A,o  A  cos  [  Am^  -  +  (j),  ] 


I  2  2 

=  A^c  ^^2  sin[AMg  +aton2(-sin  A,2cos  A)]  +  — A^c  sin[AMQ  -?Jc^o  +  ^^l] 

-A,oA  sin [ Amq  +  n^tg  -(1)^]-  A.o A  sin [ Aw^  +  A ] “ ^ ^c,o A  sin [ Am^  - +  A ] 


A  =^maxnc  COS^i^o 

Each  of  the  VC  parameters  has  a  geometric  interpretation;  the  following 
descriptions  are  based  on  my  observations.  As  with  the  circular-chief  ROEs,  and  y/^ 
are  the  amplitude  and  initial  phase  of  the  deputy's  oscillatory  out-of-plane  motion.  A 
roughly  indicates  the  scale  of  the  periodic  portion  of  the  deputy's  motion.  In  fact,  an 
upper  bound  on  the  amplitude  of  such  motion  (depending  on  phasing,  a  particular  case's 
amplitude  may  be  smaller)  is  (4/3)  Aj  -i-  A  •  A  nnd  A  nre  initial  phase  angles  for  in¬ 
plane  motion;  they  control  the  deputy's  starting  point  along  its  periodic  trajectory,  as  well 
as  the  orientation  and  "skewness"  of  that  trajectory.  A  controls  the  rate  of  drift; 
when  A  equals  0  or  ;r ,  there  is  no  drift,  and  the  drift  rate  is  a  maximum  when  A2  equals 
njl  or  3nl2.  A  scales  the  drift  rate. 

In  addition,  the  approximate  center  of  the  periodic  motion  is  (O,  A  A) ;  Eqs. 

(4.45)  reveal  why  this  is  so.  Am  (t)  -1-  nj  is  an  increasing  function  of  t ,  and  it  can  be 
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shown  (see  Appendix  G)  that  for  chief  eccentricities  less  than  about  0.3 1 ,  Am  (t )  -  n^t  is  a 
decreasing  one.  Because  of  this,  the  trigonometric  functions  in  Eqs.  (4.45)  whose 
arguments  are  Am  (t)  +  n^t  -cj)^  or  Au[t)-  nj  +  will  cycle  through  the  full  range  of 

values  between  positive  and  negative  1.  Thus,  the  first  term  may  take  on  any  value 
between  positive  and  negative  A, ,  and  the  third  term  any  value  between  positive  and 

negative  (1/3)  Aj ,  in  both  the  x-  and  y-equations.  But  because  Am  (t)  is  a  small 

oscillation  about  zero,  the  second  term  in  the  x-direction  is  always  close  to  zero,  and  the 
second  term  in  the  y-direction  is  always  close  to  .  Likewise,  this  effect  explains  why 
the  drift  term  in  the  x-direction  produces  an  oscillation  of  increasing  amplitude  about 
zero,  but  the  drift  term  in  the  y-direction  produces  a  small  oscillation  about  a  drift  rate  of 
-  (3/2)  A2  sin  • 

Bounded  Motion. 

Above,  we  eliminated  A3  from  our  parameter  set.  Let  us  examine  the  physical 
interpretation  of  A3 .  If  we  restrict  our  problem  to  bounded  motion  (equal  chief  and 
deputy  periods),  then  the  drift  terms  in  Eqs.  (4.45)  should  vanish.  In  other  words,  a  zero 
value  for  A3  represents  a  boundedness  condition  for  the  Virtual  Chief  model.  This 
condition  can  equivalently  be  expressed  in  any  of  the  following  three  ways: 

0  =  2A  +  C 
0  =  Aj  sin  ^2 

0  =  cos  [ Vc.o  -  ^c,0  ]  ( +  ^c.o )  ^0  +  sin  [ -  4^ c,o  ]  ^0 
-  sin  [vc,o  -  +  l^c.o  )  3^0  +  cos  [  Vc,o  “  ]  >'0 
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If  the  chiefs  eccentricity  is  zero,  this  reduces  to  the  HCW  boundedness  condition, 
Eq.  (3.7).  Also  note  that  Eq.  (4.49)  is  different  from  the  boundedness  condition  derived 
from  Lawden’s  equations,  Eq.  (4.3). 
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V.  Numerical  Results  and  Conclusions 


Examples 

The  Virtual  Chief  relative  satellite  motion  model  developed  in  Chapter  IV  is 
intended  to  relax  the  eccentricity  restriction  of  the  HCW  model.  That  is,  a  relative 
satellite  trajectory  predicted  by  the  Virtual  Chief  model  should  be  more  accurate  than  one 
predicted  by  HCW  for  non-zero  chief  eccentricities.  Another  way  of  stating  this 
requirement  is  that  a  Virtual  Chief  trajectory  should  be  within  a  desired  range  of  accuracy 
over  a  greater  range  of  chief  eccentricities  than  is  an  HCW  model. 

In  order  to  illustrate  the  relative  satellite  motion  predicted  by  these  models, 
example  trajectories  are  calculated  numerically  via  MATLAB  (see  Appendix  F)  and 
depicted  below.  The  trajectories  labeled  "Virtual  Chief"  are  calculated  using  Eqs.  (4.45). 
The  curve  labeled  "2BP"  follows  the  relative  position  predicted  by  unperturbed  nonlinear 
Keplerian  two-body  motion.  The  curve  labeled  "HCW"  is  calculated  by  ignoring  the 
eccentricity  of  the  chiefs  orbit  and  using  the  traditional  HCW  state  transition  matrix  (the 
first  approach  described  at  the  beginning  of  Chapter  IV).  Plots  of  the  trajectories  are 
shown  projected  into  each  c-frame  coordinate  plane  and  in  three  dimensions.  In  each 
plot,  an  "o"  marks  the  initial  position,  and  an  "x"  indicates  the  direction  of  initial  motion. 
Finally,  the  error  in  each  trajectory  (as  compared  to  the  "truth"  model,  2BP)  is  shown 
versus  time  in  each  coordinate  and  in  magnitude. 

Conditions  for  Example  I  are  shown  in  Table  1;  note  that  the  first  two  columns  are 
inputs,  and  the  third  column  is  determined  by  the  first  two  via  Eqs.  (4.45)  evaluated  at  . 
Step  size  refers  to  increments  of  chief  true  anomaly.  In  this  example,  Eigures  5  and  6,  all 
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three  models  show  negligible  drift.  It  is  clear  from  the  final  plot  in  Figure  6  that  during 


some  portions  of  the  simulated  trajectory,  the  HCW  prediction  has  a  smaller  error 
magnitude,  and  during  other  portions,  the  Virtual  Chief  prediction  has  a  smaller  error 
magnitude. 

The  error  magnitude  at  each  time  step  is  simply  the  root  sum  square  of  the  three 
error  components  at  that  time  step  (i.e.,  error  magnitude  =  x  error  +  y  error  +  z  error  ). 
It  is  also  possible  to  calculate  an  average  error  value  over  a  given  time  period  by  taking  a 
root  mean  square  of  the  error  magnitude  values  for  all  the  time  steps  during  the  given 
period  (i.e.,  for  n  time  steps,  rms  error  =  (1/n)  *  S  error  magnitude  ).  Over  the  two  chief 
orbits  plotted  for  Example  I,  the  root  mean  square  error  for  the  Virtual  Chief  trajectory  is 
15.0240191299859  meters,  and  the  root  mean  square  error  for  the  HCW  trajectory  is 
16.7326981169077  meters. 


Table  1.  Example  I  Conditions 


=8000000  (m) 

T,  =  2400  (m) 

=  -1.59999786667033  (m) 

=  0.001 

=  -K  (rad) 

=  -799.998400002667  (m) 

^c.o  =  ^  (rad) 

4  =  2400  (m) 

=  1000  (m) 

Number  of  chief  orbits:  2 

O 

II 

4  =-1.41173271214503  (m/s) 

Simulation  step  size: 

;r/360  (rad) 

^™-1000  (m) 

4  =0.00282346730660115  (m/s) 

II 

o 

5^0  =  ^  (rad) 

=  5.40274864382506e-017  (m/s) 

91 


Figure  5.  Example  I  Trajectory 
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Error  as  compared  with  two-body  solution 
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Figure  6.  Example  I  Error 
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For  some  RPO  missions,  it  is  practical  to  set  an  error  threshold  for  the  deputy; 
when  a  deputy’s  position  error  exceeds  the  threshold  distance  from  some  desired  relative 
trajectory,  a  correction  maneuver  must  be  initiated.  Figure  7  shows  a  close-in  view  of  the 
error  magnitude  during  the  first  portion  of  the  Example  I  simulation.  If  we  arbitrarily 
select  an  error  threshold  of  15  meters,  we  see  that  in  this  particular  case,  the  Virtual  Chief 
trajectory  does  not  exceed  the  threshold  until  after  approximately  10,000  seconds,  or 
about  an  orbit  and  a  half.  The  HCW  trajectory  would  exceed  this  threshold  after  about 
4000  seconds. 


Table  2.  Example  II  Conditions 


-  8000000  (m) 

> 

II 

o 

=  19.9993333633323  (m) 

=  0.005 

o 

II 

-ST 

=1999.90000416646  (m) 

(rad) 

2 

T,  =  2000  (m) 

z„=0 

Number  of  chief  orbits:  2 

O 

II 

=-0.000110280718448319  (m/s) 

Simulation  step  size: 

;r/360  (rad) 

=20  (m) 

=1.10282556487623e-006  (m/s) 

o 

II 

o 

=0.0176467162534918  (m/s) 
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Conditions  for  Example  II  are  shown  in  Table  2.  In  Figures  8  and  9,  the  Virtual 
Chief  prediction  fails  to  capture  the  drift  in  the  y-direction,  and  the  x-direction  oscillates 
in  the  opposite  direction.  In  fact,  the  trajectory  predicted  by  the  Virtual  Chief  method  is 
nearly  a  circle  confined  to  the  x-z  plane,  a  physically  unrealistic  motion.  The  Virtual 
Chief  prediction's  error  is  clearly  larger  throughout  the  simulation;  its  root  mean  square 
error  is  685.481869664407  meters.  The  root  mean  square  error  for  the  HCW  trajectory  is 
212.449649464068  meters. 


Figure  8.  Example  II  Trajectory 

Figure  10  shows  a  close-in  view  of  the  error  magnitude  plots  during  the  first 
portion  of  the  simulation.  While  in  this  case,  the  Virtual  Chief  trajectory  would  violate  a 
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15-meter  error  threshold  after  about  10  minutes,  the  HCW  trajectory  would  take  about 


twice  as  long  to  do  so. 


Figure  9.  Example  II  Error 
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Figure  10.  Example  II  Error  Magnitude 


Conditions  for  Example  III  are  shown  in  Table  3.  This  in-plane-only  example  was 


carefully  designed  so  that  the  deputy  would  follow  the  same  trajectory  as  the  chief,  as 


seen  by  the  virtual  chief.  In  other  words,  and  —  were  selected  to  be  nearly  the 

/o  dt  /o 


opposite  of  the  components  of  r^,  and  —  . 

/o  dt  /o 
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Table  3.  Example  III  Conditions 


-  8000000  (m) 

T,  =  24000  (m) 

=  -16000  (m) 

=  0.001 

7T  r  x 

^1  =  -  (radj 

2 

O 

II 

o 

^c.o  =  ^  (rad) 

4=0 

z„=0 

Number  of  chief  orbits:  2 

^2  = - (rad ) 

2 

i„=0 

Simulation  step  size: 

;r/360  (rad) 

4  =  28.2065465  (m/s) 

fo=0 

z„  =0 
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Figure  11.  Example  III  Trajectory 


Here,  as  shown  in  Figures  1 1  and  12,  the  Virtual  Chiefs  prediction  has  the  smaller 
error  throughout  the  simulation.  The  root  mean  square  error  magnitude  for  the  Virtual 
Chief  and  HCW  predictions,  respectively,  are  0.182468634831835  meters  and 
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720.129883483902  meters.  Figure  13  shows  a  close-in  view  of  the  error  magnitude  plots. 
The  HCW  trajectory  violates  a  15-meter  threshold  after  about  700  seconds,  while  in  this 
simulation,  the  Virtual  Chief  trajectory's  error  never  exceeds  half  a  meter. 


Error  as  compared  with  two-body  solution 
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Figure  12.  Example  III  Error 


Figure  13.  Example  III  Error  Magnitude 
Error  Analysis:  Phase  1 

We  have  already  established  that  every  dynamics  model  represents  a  trade-off 
between  simplicity  and  accuracy.  Since  the  Virtual  Chief  method  is  more  complex  than 
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the  HCW  model,  it  is  worth  examining  what  we  have  gained  in  accuracy.  Examples  such 


as  those  in  the  previous  section  begin  to  suggest  an  answer. 

Numerical  simulations  reveal  that  in  some  cases,  the  HCW  model  predicts  the 
relative  motion  with  better  accuracy  than  the  Virtual  Chief  method.  To  illustrate  how  the 
comparative  accuracy  changes  by  case,  I  conducted  a  series  of  simulations  using  the  same 
Virtual  Chief  parameters  as  Example  II  above.  In  each  case,  the  chiefs  semi-major  axis 
was  16,778.137  km.  The  chiefs  eccentricity  and  mean  anomaly  at  epoch  were  varied. 


Table  4.  Model  Accuracy,  Series  1 


n 

4 

n 

2 

In 

4 

Ik 

8 

n 

9k 

8 

5k 

4 

3;r 

2 

13;r 

8 

Ik 

4 

15;r 

8 

^c=0 

tie 

tie 

tie 

tie 

tie 

tie 

tie 

tie 

tie 

tie 

tie 

Tie 

0.00010 

H 

H 

H 

H 

H 

V 

V 

V 

V 

V 

V 

V 

0.00015 

H 

H 

H 

H 

H 

V 

V 

V 

H 

H 

V 

V 

0.00020 

H 

H 

H 

H 

H 

V 

V 

H 

H 

H 

H 

V 

0.00100 

H 

H 

H 

H 

H 

V 

H 

H 

H 

H 

H 

H 

0.00500 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

0.01000 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

Table  4  summarizes  the  results,  showing  which  model  was  more  accurate,  as 
measured  by  root  mean  square  error.  Each  column  is  a  different  value  of  „  and  each 

row  a  different  value  of  .  An  "H"  indicates  that,  for  that  case,  the  HCW  prediction  had 

the  lower  root  mean  square  error  over  the  course  of  the  simulation.  A  "V"  indicates  that 
the  Virtual  Chief  prediction  had  the  lower  root  mean  square  error.  As  previously  noted. 


98 


when  is  zero,  the  two  models  are  identieal.  Inspeetion  of  Table  4  seems  to  indieate 


that,  for  these  examples,  the  Virtual  Chief  method  performed  worse  as  chief  eccentricity 
increased,  and  that  it  only  performed  better  than  HCW  during  the  half-orbit  from  chief 
apogee  to  perigee. 

In  essence,  I  had  been  holding  the  relative  motion  problem  constant  and  varying 
the  chief  orbit  conditions  to  find  some  region  where  the  Virtual  Chief  model  represents  a 
good  choice.  Since  Table  4  does  not  suggest  such  a  region,  I  attempted  another  series  of 
simulations  holding  the  chief  orbit  conditions  constant  and  varying  the  in-plane  Virtual 
Chief  parameters. 

This  second  series  of  simulations  used  a  chief  eccentricity  of  0.00015  and  a  chief 
mean  anomaly  at  epoch  of  ,  corresponding  to  a  case  in  Table  4  for  which  the  Virtual 
Chief  model  performed  better.  The  results  are  presented  in  Table  5. 


Table  5.  Model  Accuracy,  Series  2 


A  (m) 

A  (rad) 

A  (m) 

A  (rad) 

“Winner” 

0 

0 

0 

0 

Tie 

0 

0 

2 

0 

H 

2 

0 

0 

0 

H 

2 

0 

2 

0 

H 

0 

0 

2 

n/ 

74 

V 

2 

n/ 

74 

0 

0 

V 

2 

0 

2 

n/ 

74 

H 

2 

n/ 

74 

2 

0 

V 

2 

n/ 

74 

2 

n/ 

74 

H 
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Thus,  for  the  cases  in  Table  5,  the  Virtual  Chief  model  did  not  consistently  out¬ 
perform  the  HCW  model.  The  results  of  Tables  4  and  5  imply  that  searching  either  the 
parameter  space  of  the  chiefs  orbit  elements  or  the  parameter  space  of  the  relative  initial 
conditions  (here,  in  the  form  of  the  Virtual  Chief  parameters)  will  not  find  a  region  of 
predictably  good  performance  for  the  Virtual  Chief  model. 

We  can  make  a  few  statements  about  the  modeling  error  in  this  study.  The 
models  compared  (the  two-body  truth  model,  HCW,  and  Virtual  Chief)  have  three 
possible  error  sources:  unmodeled  forces,  eccentricity  effects,  and  nonlinear  effects. 
They  correspond  to  the  three  assumptions  underlying  the  equations  of  motion.  First,  the 
two-body  assumption  incurs  errors  from  unmodeled  differential  forces  (J2  and  higher- 
order  gravity  terms,  differential  drag,  etc.).  This  error  is  acceptable  over  short  time 
scales,  and  is  common  to  HCW,  Virtual  Chief,  and  the  two-body  truth  model.  Therefore, 
this  error  source  has  no  effect  in  the  current  study. 

Second,  the  circular  chief  assumption  inherent  to  the  HCW  model  incurs  error 
because  the  chiefs  distance  to  the  Earth  and  angular  rate  are  not  constant.  This  is  the 
error  which  the  Virtual  Chief  method  was  designed  to  remove  or  lessen.  The  two-body 
truth  model  does  not  incur  this  error.  In  the  HCW  model,  this  error  can  be  approximated 
as  an  increasing  function  of  actual  chief  eccentricity,  (^c )  • 

Third,  linearization  of  the  gravity  field  incurs  errors  because  the  higher-order 
relative  position  terms  are  neglected  ( ,  xy ,  etc.).  The  two-body  truth  model  does  not 
incur  this  error.  In  the  HCW  model,  this  error  should  tend  to  increase  with  the 
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displacement  between  deputy  and  chief.  Thus,  if  is  an  increasing  function  of 


distance  and  p  is  ,  we  can  approximate  the  HCW  error  as 

4cw  ~  >1  [fecc  (^’c ) ’  fdisp  (P))  (5-1) 

The  Virtual  Chief  model  is  not  immune  to  linearization  error;  in  a  sense,  we  are 
trying  to  fold  the  eccentricity  error  into  the  displacement  error.  The  Virtual  Chief 
displacement  error  has  two  contributions,  because  we  linearized  both  the  displacement 
from  the  actual  chief  to  the  virtual  chief,  ,  and  the  displacement  from  the  deputy  to  the 

virtual  chief,  .  Obviously,  the  error  component  dependent  on  would  vary  by 
individual  case. 

Thus,  we  can  approximate  the  Virtual  Chief  error  as 

Svc  ~  f^i  [fdisp  (Pc  )  ’  fdisp  (Pd  ))  (5-2) 

Increasing  tends  to  increase  ,  but  also  tends  to  increase  p^ ,  and  therefore 
has  an  effect  on  .  This  effect  would  be  consistent  with  the  behavior  seen  in  Table  4. 

It  may  also  prevent  the  Virtual  Chief  method  from  being  used  to  relax  the  HCW  model's 
eccentricity  restriction,  at  least  in  the  general  case. 

Also,  since  the  chief  is  in  an  orbit  coplanar  with  that  of  the  virtual  chief,  and  since 
no  along-track  offset  is  present,  the  separation  between  the  chief  and  the  virtual  chief  is 
due  entirely  to  the  eccentricity  difference.  In  fact,  from  the  HCW  in-plane  solution  to  the 
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circular  chief  problem,  we  know  that  the  chiefs  motion  relative  to  the  virtual  chief 
describes  a  2x1  ellipse  with  maximum  separation  occurring  in  the  along-track  direction. 
The  value  of  this  maximum  displacement  is  found  from  Eq.  (3.32)  to  be  2a^e^ .  This 
tends  to  be  a  very  large  quantity,  even  for  very  low  eccentricities  and  semi-major  axes.  If 
the  deputy  is  far  away  from  the  chief,  so  that  p  approaches  the  value  of  'la^e^ ,  then 
is  likely  also  very  large.  This  would  only  be  untrue  if  the  trajectory  is  carefully  designed 
so  that  the  deputy  remains  as  close  to  the  virtual  chief  as  does  the  actual  chief,  as  in 
Example  in. 

Eor  example,  the  conditions  used  to  generate  Table  4  gave  the  deputy  a  large 
positive  along-track  displacement  from  the  chief;  in  other  words,  the  deputy  is  leading  the 
chief.  During  the  half-orbit  from  perigee  to  apogee,  the  chief  is  itself  leading  the  virtual 
chief,  leading  to  a  large  displacement  between  the  deputy  and  the  virtual  chief.  However, 
during  the  half-orbit  from  apogee  to  perigee,  the  virtual  chief  is  leading  the  chief,  placing 
the  virtual  chief  in  the  vicinity  of  the  deputy.  This  may  explain  why,  for  the  example  of 
Table  4,  the  Virtual  Chief  prediction  could  only  be  more  accurate  during  the  second  half¬ 
orbit. 

Error  Analysis;  Phase  2 

If  there  were  some  way  to  define  a  new  parameter  space  (other  than  chief  orbit 
elements  and  relative  initial  conditions)  in  which  we  could  easily  specify  that  the  deputy 
must  remain  relatively  close  to  the  virtual  chief,  then  it  may  be  possible  to  discover  a 
restricted  region  of  good  accuracy  for  the  Virtual  Chief  model.  As  a  test,  let  us  use  the 
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deputy’s  ROEs,  as  seen  from  the  virtual  chief;  that  is,  the  ROEs  we  would  use  to  describe 
the  HCW  motion  of  the  deputy  in  the  o-frame.  These  deputy  o-frame  ROEs  represent  a 
combination  of  the  information  in  both  the  chiefs  orbit  elements  and  the  deputy’s 
relative  initial  conditions. 

Note  that  the  chief’s  o-frame  ROEs  are  always  zero  except  for  and  p .  Note 
also  that  the  deputy’s  is  the  same  for  both  the  o-frame  HCW  model  and  for  the  c- 
frame  Virtual  Chief  model. 

Eor  the  sake  of  simplicity,  let  us  hold  constant  the  chiefs  orbit  using  the  same 
conditions  as  Example  HI  above  (Table  3)  and  examine  bounded,  in-plane-only  relative 
motion.  Let  us  further  restrict  the  deputy  to  o-frame  motion  centered  about  the  virtual 
chief,  so  that  the  only  ROEs  which  can  be  non-zero  are  and  P .  As  shown  in  Eigure 
14,  these  restrictions  are  equivalent  to  varying  only  the  deputy’s  c-frame  value,  while 
constraining  to  equal  -2n^X(, .  In  the  figure,  the  various  dots  along  the  x-axis  are 

some  of  the  sampled  deputy  initial  positions.  Note  that,  under  these  restrictions,  the  chief 
is  always  at  apogee  at  epoch.  Likewise,  the  deputy  is  at  apogee  or  perigee  at  epoch. 

By  varying  the  deputy’s  Xg  under  these  conditions,  we  are  effectively  varying  its 

.  The  ratio  of  the  deputy’s  o-frame  to  the  chiefs  o-frame  gives  an  indication  of 

the  relative  closeness  of  the  deputy  to  the  virtual  chief,  which  we  expect  to  improve  the 
Virtual  Chief  modeling  error.  In  Table  6  below,  the  index  column  simply  represents  this 
ratio. 
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T 


Figure  14.  Phase  2  Geometry 


Tables  6  give  the  simulation  results  for  this  phase  of  the  study.  All  values  are  in 
meters,  except  for  the  index.  The  tables  are  organized  by  increasing  index  value.  Note 
that  when  the  index  is  zero,  the  deputy  is  co-located  with  the  virtual  chief.  The  other 
table  entries  occur  in  pairs;  this  is  because  the  index  will  have  the  same  value  regardless 
of  whether  the  deputy  is  at  perigee  or  apogee. 

Note  that  the  Virtual  Chief  model  performed  with  accuracy  equal  to  or  better  than 
HCW  for  most  of  the  cases  in  Table  6.  These  correspond  to  deputy  o-frame  values 
both  less  than  and  greater  than  the  chiefs  o-frame  ;  in  other  words,  for  deputy 
eccentricities  both  less  than  and  greater  than  the  chiefs  eccentricity.  Since  this  one¬ 
dimensional  search  indicates  a  region  of  predictable  accuracy  for  the  Virtual  Chief  model. 
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we  can  begin  to  look  for  an  analytical  understanding  of  the  modeling  error  by  capturing 
the  actual  root  mean  square  error  values  (the  last  two  columns  of  Table  6). 


Table  6.  Phase  2  Results 


^0 

Chief 

Deputy 

Index 

VC  rms  error 

HCW  rms  error 

-8000 

16000 

0 

0 

90.76991438 

450.51339369 

-6000 

16000 

4000 

0.25 

85.00250365 

354.74757485 

-10000 

16000 

4000 

0.25 

85.22566297 

535.01964956 

-12000 

16000 

8000 

0.5 

68.35298262 

608.25772993 

-4000 

16000 

8000 

0.5 

67.94038372 

247.73234487 

-14000 

16000 

12000 

0.75 

40. 1 36294 1 1 

670.22102227 

-2000 

16000 

12000 

0.75 

39.60054922 

129.47905476 

0 

16000 

16000 

I 

0 

0 

-16000 

16000 

16000 

I 

0.182468635 

720.12988348 

2000 

16000 

20000 

1.25 

50.84427082 

I40.69I77273 

-18000 

16000 

20000 

1.25 

50.42206374 

760.31085429 

-20000 

16000 

24000 

1.5 

II2.7794I75 

788.44127984 

4000 

16000 

24000 

1.5 

II2.9I528I9 

292.58261425 

6000 

16000 

28000 

1.75 

I86.I960657 

455.65838792 

-22000 

16000 

28000 

1.75 

186.5534288 

805.3078I9I8 

-24000 

16000 

32000 

2 

271.7596827 

810.93090792 

8000 

16000 

32000 

2 

270.6696696 

629.90456028 

16000 

16000 

48000 

3 

720.1537784 

1438.2923102 

-32000 

16000 

48000 

3 

727.2472372 

722.05843779 

Figures  15-18  illustrate  the  c-frame  x-y  plane  trajectory  predictions  for  four  cases 
from  Table  6.  It  is  clear  that  HCW  is  always  predicting  bounded  motion  (repeating 
trajectories).  The  accuracy  of  the  Virtual  Chief  prediction  seems  to  depend,  for  these 
cases,  on  how  well  it  captures  the  drift  in  the  truth  model.  Figure  18  corresponds  to  the 
only  case  in  Table  6  for  which  HCW  showed  superior  accuracy. 
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Relative  trajectory 


x(m) 


Figure  15.  Index  =  0 


xio*  Relative  trajectory 


Figure  16.  Index  =  1 
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xio*  Relative  trajectory 


Figure  17.  Index  =  3  (xO  =  16  km) 


xio^  Relative  trajectory 


Figure  18.  Index  =  3  (xO  =  -32  km) 
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We  can  also  plot  the  root  mean  square  error  for  each  model  against  the  deputy’s  c- 
frame  Xq  ,  as  in  Figure  19.  Also  depicted  on  the  figure  are  analytic  functions  heuristically 

matched  to  the  root  mean  square  error  data.  For  the  cases  in  Phase  2,  both  the  HCW  error 
and  the  Virtual  Chief  error  functions  behave  as  a  quadratic  form: 


'HCW 


'VC 
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Figure  19.  Phase  2  Error  versus  xO 


Figure  20  shows  a  similar  plot  of  the  Virtual  Chief  error  as  a  function  of  the  Table 
6  index  value.  Again,  a  heuristically  derived  quadratic  function  matches  the  error  data: 
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Syf.  =  90  *  abs  index^  - 


Figure  20.  Phase  2  VC  Error  versus  Index 


One  result  from  Phase  2  that  is  perhaps  surprising  is  that  dy^  is  not  an  increasing 
function  of  Pjj  (tp) .  Rather,  the  error  is  lowest  when  the  ratio  is  closest  to  unity. 

What  we  have  achieved,  though,  is  a  region  of  predictable  accuracy  for  the  Virtual  Chief 
model.  It  is  clear  from  Figure  19  that,  subject  to  the  Phase  2  restrictions,  for  any  value  of 
Xq  greater  than  about  -32,000  meters,  the  Virtual  Chief  root  mean  square  error  should  be 
less  than  the  HCW  root  mean  square  error. 


Conclusions  and  Recommendations 

Relative  Orbit  Elements  constitute  a  useful  realization  of  the  HCW  model  for  the 
circular  chief  problem.  As  we  saw  in  Chapter  HI,  they  provide  significant  geometrical 
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insight  into  the  relative  motion.  ROEs  merit  continued  development  to  seek  relative 
navigation  techniques,  autonomous  maneuver  planning  algorithms,  and  other  RPO 
applications. 

The  latter  part  of  this  study  has  developed  an  approximate  relative  satellite  motion 
model  for  chief  and  deputy  with  slightly  eccentric  orbits.  This  model  has  a  simple 
parameterization  designed  for  a  variety  of  applications,  including  autonomous  RPO.  The 
model  achieves  this  simplicity  at  the  risk  of  lower  accuracy;  in  some  cases,  using  the 
HCW  model  may  be  a  better  choice.  In  fact,  in  some  cases,  the  new  model  tends  to  break 
down  with  increasing  eccentricity. 

The  overarching  goal  was  to  find  a  simple  parameterized  model  that  provides 
geometrical  insight  and  operational  efficacy  while  predicting  relative  motion  accurately  at 
chief  eccentricities  higher  than  can  the  HCW  model.  Under  the  restrictions  for  the  Phase 
2  error  analysis,  this  has  been  achieved  for  a  predictable  region.  Further  study  is  clearly 
indicated  to  expand  this  region  of  predictable  accuracy.  Specifically,  we  should  search 
the  deputy’s  o-frame  ROE  space  under  different  cases,  varying  multiple  parameters  and 
characterizing  the  behavior  of  the  root  mean  square  error  for  both  the  Virtual  Chief  and 
the  HCW  models. 

However,  the  broader  goal  of  finding  a  simple  model  that  is  generally  more 
accurate  than  HCW  for  elliptical  orbits  will  have  to  be  reached  via  a  different  approach. 
Any  of  the  possible  approaches  listed  at  the  beginning  of  Chapter  IV  remain  open.  With 
further  study,  these  other  options  may  prove  a  viable  answer. 
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For  example,  the  quadric  surface  of  relative  motion  found  by  Jiang,  et  al.  [45], 
might  be  adapted  for  drifting  cases,  so  that  the  quadric  surface’s  parameters  are  time- 
varying.  Just  as  HCW  predicts  a  drifting  instantaneous  ellipse,  Lawden’s  equations  may 
predict  a  drifting,  oscillating  instantaneous  quadric  surface. 

Also,  it  may  be  possible  to  find  geometrically  meaningful  parameters  in  a  first- 
order,  time-explicit  approximate  model,  such  as  one  derived  from  the  first-order  state 
transition  matrix  of  Melton  [66] . 

Alternatively,  a  linearized  relative  satellite  motion  model  expressed  in  coordinates 
of  a  non-rotating  rectangular  reference  frame,  fixed  to  the  chief  and  always  parallel  to  the 
chiefs  perifocal  frame,  may  prove  to  be  simpler  than  existing  LVLH-frame  models. 

We  have  already  examined  the  fifth  approach  listed  in  Chapter  IV,  the  Virtual 
Chief  Method,  in  its  linearized  form.  However,  if  a  similar  method  were  used  to  develop 
a  quadratic  model  (second-order  in  the  relative  distances),  the  accuracy  problems 
discovered  in  Chapter  V  would  likely  improve.  This  would  be  very  close  to  the  approach 
of  Fasano  and  D’Errico  [30].  However,  this  accuracy  improvement  would  come  at  the 
expense  of  simplicity. 
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Appendix  A:  Identities 


Common  Trigonometric  Identities 


For  convenience,  the  trigonometric  identities  used  in  this  thesis  are  listed  here: 


sin^a+cos^a  =  l 


(A.l) 


sin(a  +  P)  =  cos  a  sin  /? + sin  a  cos  P 


(A.2) 


sin  a  =  sin(;r-a) 


(A.3) 


cos  a  =  -  cos 


(;r-a) 


(A.4) 


i(-a)  =  - 


(A.5) 


=  cos« 


(A.6) 


cosacosy^  =  ^|cos(a-y0)  +  cos(a  +  y0)| 


(A.7) 


sin«siny5  =  — |cos(« -y0)-cos(« +  y0)| 


(A.8) 


cosasin  P  =  ^|sin(a  +  y0)-sin(a-y5)| 


(A.9) 
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(A.  10) 


sin«cosy0  =  ^|sin(«  +  y5)  +  sin(«-y5)| 


cos(a  +  P)  =  cos  a  cos  /?  -  sin  a  sin  /? 


(A.ll) 


cos(;r  +  a)  =  -cos  (a) 


(A.12) 


Identities  for  the  atan2  Function 

Similar  identities  for  atan2,  the  quadrant- specific  inverse  tangent  function,  are  less 
commonly  found.  The  identities  used  in  this  thesis  are  listed  below.  They  can  often  be 
understood  graphically  from  right  triangles  and  a  unit  circle,  as  shown  in  Figure  21. 

atan2[-p,a)  =  -atan2[p,a)  (A.  13) 

atan2[p,-a^  =  TT -atan2{^P,a^  (A. 14) 

atan2[-/3,-a)  =  TT  +  atan2[/],a)  (A. 15) 
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Appendix  B:  Harmonic  Addition  Theorem 

To  deal  with  periodic  terms  in  the  trajectory  equations,  I  used  the  principle  stated 
in  Eqs.  (3.9)  through  (3.11)  and  repeated  below,  sometimes  referred  to  as  the  Harmonic 
Addition  Theorem  [101].  A  periodic  function  defined  as 

f{6)  =  Asm6+B(x>^6 

can  be  equivalently  written  both  as 

f{6)  =  +B^  sm(d +  atan2(B,A)) 


and  as 


f(0)  =  -yjA^  +B^  cos(0  +  atan2(A,-B)) 

where  atan2  is  the  quadrant- specific  inverse  tangent  function. 

The  first  step  in  deriving  Eqs.  (3.10)  and  (3.11)  from  Eq.  (3.9)  is  to  find  the 
amplitude  of  the  oscillation  described  by  Eq.  (3.9).  Since  Eq.  (3.9)  is  written  without  an 
offset  term,  the  oscillation  is  centered  about  zero,  and  the  peak  values  of  f(0)  will  equal 
the  amplitude.  Each  peak  value  is  a  local  maximum  where  the  derivative  of  f(0)  equals 
zero: 
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f{0)  =  0  =  Acos  O-BsinO 


Rearranging, 


Acosd  =  Bsind 


Squaring  both  sides, 


A"cos"6»  =  fi"sin"6» 


(B.l) 


Substituting  the  trigonometric  identity  in  Eq.  (A.l)  into  the  right-hand  side,  Eq. 
(B.l)  becomes 


A^  cos^  6  =  B^  -B^  cos^  0 


Solving  for  the  constant  term. 


=(AVb")cos"6' 


Taking  the  square  root. 


B  =  ±Va"  +  B"  cosB 


Rearranging, 
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±B 


(B.2) 


cos  6  =  , - 

Substituting  the  identity  from  Eq.  (A.l)  into  the  other  side  of  Eq.  (B.l)  and  then 
repeating  the  same  steps  as  above  yields  the  result 

(B.3) 

Eqs.  (B.2)  and  (B.3)  define  the  extrema  of  f{6) .  If  we  choose  the  positive  value 
from  each,  choose  6^^  as  the  corresponding  value  of  6 ,  and  substitute  into  Eq.  (3.9), 
we  can  then  write 


fiO  J  =  A( 

'  m«Y  '  ' 


)  +  5(- 


B 


0 


where  )  represents  the  peak  value,  the  amplitude  we  have  been  seeking. 

Simplifying, 


A^  +  B^ 
+  B^ 


Now  that  we  know  the  amplitude  of  f{6) ,  we  can  convert  Eq.  (3.9)  into  a  single 
sinusoidal  term  by  factoring  out  the  amplitude: 
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f{e)  =  ^A^  +  B^ 


A  ■  n  ^  n\ 

rSin6^+  , - ^=cos6^ 


1  UiXX  I  i 


(B.4) 


The  coefficients  inside  the  braces  suggest  a  right  triangle  with  legs  of  length  A 
and  B.  Its  hypotenuse  would  have  length  +  B^  .  One  of  this  triangle's  angles  would 

have  a  cosine  equal  to  ,  ,  a  sine  equal  to  ,  ,  and  a  tangent  equal  to  -  . 

A 


If  we  identify  this  angle  as  atan2{B,A) ,  then  we  can  write  Eq.  (B.4)  as 

f{0)  =  ^Ia^  +  B^  {cos(atan2(B,  A)) sin 6*  +  sin(aton2(B,  A))cos 6*}  (B.5) 

Using  the  trigonometric  identity  from  Eq.  (A. 2),  this  becomes 

f(d)  =  ^Ja^  +B^  sin(0  +  atan2(B,  A)) 

which  is  Eq.  (3.10). 

Next,  we  can  substitute  Eq.  (A.  16)  into  Eq.  (3.10),  yielding 


f(0)  =  ^/A"+B"  sin(^  +  atan2(A,  -B)  -  ^) 


(B.6) 


If  we  apply  the  trigonometric  identities  from  Eqs.  (A. 5)  and  (A. 6),  Eq.  (B.6)  becomes 
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f(0)=-4^7¥  cos  (^  +  atan2(A,  -B)) 


which  is  Eq.  (3.11). 

Both  of  these  forms,  Eqs.  (3.10)  and  (3.11),  are  periodic  functions  whose 


amplitude  (\Ia^  +  B^  ),  initial  phase  (atan2(B,A)  or  atan2(A,—B)),  and  time-varying 
phase  parameter  ( )  are  all  readily  apparent. 
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Appendix  C:  An  Ellipse  Oriented  Arbitrarily  in  Three  Dimensions 

In  general,  to  define  a  curve  or  trajectory  in  3-D  space  requires  at  least  two 
equations.  If  the  curve  is  an  ellipse,  then  we  can  define  a  principal- axis  reference  frame 
with  unit  vectors  u  ,v ,  and  w  .  In  this  frame,  u  and  v  are  aligned  with  the  ellipse's  major 
and  minor  axes  (although  it  is  not  yet  necessary  to  specify  which  is  which),  and  w  is 
normal  to  the  plane  containing  the  ellipse. 

Let  position  with  respect  to  the  ellipse  center  be  defined  by  a  vector  with 
principal- axis  frame  coordinates  x',  y',  and  z': 

To/  =  x'u  +  y'v  +  z'w 

Then  the  two  equations  necessary  to  define  the  ellipse  are 

z'  =  0 


Or,  writing  the  system  parametrically, 

x  =  a  cos  6 
y'  =  b^m  6 

where  a  and  b  are  the  semi-major  and  semi-minor  axis  lengths  and  6^  is  a  phase 
parameter. 
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Now  let  us  allow  the  center  of  the  ellipse  to  translate  from  the  origin  of  the 


principal- axis  frame  to  some  point  C  whose  principal- axis  coordinates  are 


C  = 


c. 


The  form  of  the  elliptical  trajectory  is  therefore 


■x' 

a  cos  6  -1-  C^. 

y' 

= 

bsmO  +  C^, 

z'  \ 

_  C,  _ 

In  order  to  express  this  trajectory  in  arbitrary  coordinates  x,  y,  and  z,  we  must 
construct  a  rotation  matrix  by  writing  the  unit  basis  vectors  of  the  principal-axis  frame 
(u  ,v ,  and  vv )  in  terms  of  their  x-,  y-,  and  z-coordinates.  The  system  will  then  be 


"^1 

1  ^  1 

x' 

t 

1 

1 _  . 

=  [m  1  V 

1 

\ 

1 _ z _ 

Let  X-,  y-,  and  z-coordinate  values  be  specified  by  subscripts.  Then  the  parametric 
equations  of  an  elliptical  trajectory  in  an  arbitrary  3-D  reference  frame  are 


x(t)  =  au^  cos  6  +  bv^  sin  0  +  C^ 
y(t)  =  aUy  cos  6  +  bv^  sin  O  +  C^ 
z(t)  =  au^  cos  6  +  bv^  sin  0  +  C^ 


121 


Any  trajectory  equations  which  satisfy  this  form  must  therefore  describe  an 

ellipse. 
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Appendix  D:  Two-body  Orbit  Properties 

The  following  properties  of  Earth  satellite  orbits  under  the  two-body  assumption 
can  be  found  in  most  astrodynamics  texts.  They  are  listed  here  for  convenience,  as  well 
as  to  express  them  in  notation  consistent  with  this  thesis. 

The  period  P  of  an  orbit  is  related  to  the  Earth's  gravitational  parameter  fj.  and 
the  orbit's  semi-major  axis  a  as: 


(D.l) 


A  satellite's  orbit  radius,  or  its  distance  from  the  Earth's  center,  at  time  t  is  r  (t) . 
It  can  be  found  from  the  orbit's  semi-major  axis,  the  orbit's  eccentricity  e ,  and  the 
satellite's  true  anomaly  v{t)  as: 


r 


(<) 


l-l-ecosv(t) 


(D.2) 


A  satellite's  mean  anomaly  M  (t)  is  defined  in  terms  of  its  mean  motion  n  and 
time  of  perigee  passage  as: 


(D.3) 
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Kepler's  equation:  a  satellite's  mean  anomaly  is  related  to  its  eccentric  anomaly 


E{t)  as: 


M  (t)  =  £■(?)- e  sin  £■(?) 


(D.4) 


Eccentric  anomaly  is  related  to  true  anomaly  as: 


v(t)  =  latanl 


sin 


E{t) 

2 


(D.5) 


A  satellite’s  orbit  radius  at  perigee,  ,  is  found  as: 


r^=a{l-e) 


(D.6) 


The  semilatus  rectum  p  of  a  satellite’s  orbit  is  found  as: 


p  =  a(l-e^^ 


(D.7) 


The  angular  rate  v(t)  of  a  satellite’s  orbit  is  found  as: 


=  2/""  ,.2  +  (D.8) 

a  \\-e  1 

=  n(l-e^^  (l  +  ecos[v(t)]j 
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The  angular  momentum  h  of  a  satellite’s  orbit  is  found  as: 


(D.9) 


It  is  also  true  that 


h=  r(t)  ^v{t) 


(D.IO) 


Combining  Eqs.  (D.9)  and  (D.IO),  and  then  using  Eqs.  (D.2),  (D.7),  and  (D.8), 


na 


(D.ll) 


Kepler’s  equation  (Eq.  (D.4))  may  be  replaced  by  a  Eourier-Bessel  series 
expansion  of  the  true  anomaly  [6],  written  as: 


00  1 

+2y-< 


00 

Y 

/  1  / 

1 

<N 

1 

1 

1 _ 

|A:+« 

/l=-00 

y  sin  (t)) 


(D.12) 


A  satellite's  mean  motion  is 


2n 
n  =  — 
P 
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Appendix  E:  Elements  of  the  Virtual  Chief  State  Transition  Matrix 

The  sixteen  in-plane  elements  of  the  state  transition  matrix  in  Eq.  (4.36)  are  as 
follows: 


=  cos  \vc{t)-M^  (t)]  [-3cos  [n^  {t  -  to )]  +  4]  cos  [v^.o  “  4^c,o] 

[6  sin 1, )]  -  6^^  (t  -  to )]  cos  [v^.o  “  4^c,o  ] 


-l-sm 


■(^C,0  ^c) 


[+sin[vc.o-Mc.o] 

[  sin[nc(t-to)]  . 


cos 


_^C,0  4^C,o] 


+  - 


2cos[nc(t-to)]  2 


-I - Icos 

n 


-l-sm 


m\vc{t)-Mc{t)\\ 


2  cos  [^^(t- to)]  2 


4sin[nc(t-to)] 


\_^C,0  4^C,o] 

\_^C,0  “4^C,o] 

3t  j"  cos  []  Vp  0  “  CO  ] 


>sm 
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+  sin 


(0  ^ c  (0_  1 


®  12  =  -cos  [vc  (0  - (0]  [-3 cos  [n^  (?  - )]  +  4] sin  [v^  ^  ] 

_6  sin  \n^{t- 1, )]  -  [t  -  ?o )]  sin  [v^.o  “ ^c,o  ] 

+  cos[vc,o-^c,o] 


■(^C,0  ^c) 


cos 


[^c(0-^c(0] 


-cos 


ycfi  ^c,o_ 


-l-sin 


\yc  (0  ^ c  (0_  1 


tic 

-  +  — 
Tie 

2cos[nc(?-?o)] 

.21 

tie 

«cj 

4sin[nc(?-?o)] 

Tie 

>sin[ 

-3(?-?o) 

cos 


_^C,0  ^c,o_ 


^c,o_ 


sin[nc(?-?o)] 


Oi4=cos[vc(0-^c(0] 


cos 


[^C,0  -^C.o] 

“c 

2cos[nc(?-?o)]  2 


+  sin 


in[vc(0-^c(0] 


2cos[n(.(?-?o)]  2 


cos 


+ 


4sin[nc(/‘-?o)] 


-3(?-?o) 


sin[vc.o-^c.o] 

[^c,o  c,o~\ 

[^c,o  ~^c,o~\ 


sin 


[  sin[nc(?-?o)] 


Oi5=cos[vc(0-^c(0] 


sin 


[^c,o  -^c.o] 


+ 


+  sin 


[^c  (0  ^ c  (0_  I 


+ 


Tic 

He 

2cos[n(,  (?-?o) 

]  2], 

Tic  _ 

4sin[n^  (^-^o) 

l-3(?- 

cos 


sin 


\y  C,Q  -^C.o] 

[^c,o  — -^c.o] 

[^C,0  Cfi~\ 


cos 
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O21  =-sin[vc(0-^c(0_ 

[-3 cos  [hc  (?  -  ?o )]  +  4]  cos  [ Vc.o  -Mc,o_ 

+  cos[vc(0-^c(0]' 

[6  sin  [hc  (?  -  ?o )]  -  6nc  (?  -  ?o 

+  sin[vc,o-Mc,o] 

cosj^v^Q  4fc,o] 

-sin 


+ 


+  COS 


['^C  (0  (0_  1 


+  COS 


[vc{t)-M^{t)\\ 


- = - ^sin 

L'^C.O  - 

He 

2cos[np  (^“^0) 

tie 

tie 

2cos[nc(?-?o)] 

2 

tie 

«cj 

4sm[ne{t-to)] 

tic 

cos[v 

_-3(?-?o) 

- 

?o)]  +  4]sin[vco- 

Mc,o] 

“^0)]  “^0) 

sin  [v, 

cos 


sin 


\yc,Q  -^C.o] 
[^c,o 


’['^C.O  ^c,o_ 


+  COS 


■('^C,0  ^c) 


\ycfi  -^c.o] 


sin 


['^c  (0  '^c  (0_  1 


cos 


_^C,0  ^c,o_ 


+ 


-cos 


K(0-^c(0] 


tic 

F  — 
ric 

2cos[nc(?-?o)] 

2 

tic 

tic 

4sm[ne{t-to)] 

tic 

sin 

1 - 

1 

U) 

1 

0 

1 

cos 


‘  _^C,0  ^ c,o_ 
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O24=-sin[vc(0-^c(0] 


cos 


+ 


[^c,o  -^c.o] 

“c 

2cos[nc(^-^o)]  2 


-  +  - 


+  COS 


[v,{t)-M,{t)\\ 


2cos[nc  (?-?o)]  2 


cos 


^[^c,o  -^c.o] 

[^c,o  ~^c,o\ 


+ 


-3(?-?o) 


sin 


_^C,0  -^C.o] 


O25  =-sin[vc(0-^c(0_ 


sin 


+ 


“c 

2cos[np  (?-?„)]  2 


cos 


+  COS 


2cos[n^  (?-?o)]  2 


sin 


[^c,o  -^c.o] 

[^c,o  c,q\ 


+ 


4sin[nc(?-?o)] 


cos 


_^c,o  -^c.o] 
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®4,  =  COS  [v,  (0  -  A?,  (0]  [3«e  sin  [«c  )]]  cos  [v^  „  -  „  ] 

+  sin  [v^  (f )  -  (f )]  [6n,  cos  [«,(?-  f„ )]  -  6nJ  cos [v^^„  -  ] 


+  (’^c,o-«c)i 


cos 


K(0-^c(0]' 


-cos 


+  sin[v^(f) -Me  (0] 


K(^-^o)]sin[ve,o-Me  J 
+2 sin  [n^{t-  )] cos  [ve.o  - 
\ 2 sin  [n^{t-  io )]  sin  [ve.o  -  ] 

-  [4cos  [n^{t- )]  -  3]  cos  [ve.o  “  A^e.o  ] J 


-('^c(0-«c)cOs[Ve(0-A?e(0] 


r[6sin[ 


-('^c(0-«c)('^c,o-«c)j 


-sin[Ve(0-A?e(0]l 


+  COs[Ve(f)-A?e(0]l 


)]  +  4]cos[ve„-M, 

c.o  ] 

)]-6ne(f-f„)]cos 

1 

Ko- 

0  J 

sin[«c(f-fo)]  .  r 

.^c.o  - 

Sin 

«c 

2cos[ne(f-f„), 

1+- 

”c 

”c 

2cos[ne(f-fo)] 

2 

”c 

”c_ 

4sin[ne(f-fJ] 

”c 

cos[i 

COS 


inKo-A?c,o] 
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®  42  =  -  cos  [v,  (f )  -  (f )]  [3n,  sin  [«,(?-?„)]]  sin  [v,  „  -  „  ] 

-  sin  [^c  (0  -  (0]  [K  cos  [n,  (f  -  )]  -  6n,  ]  sin  ] 


—  (v  )  1 

\  c,n  c  / 


COS 


cos 


+  sin 


K(0-^c(0] 

inK(0-^c(0]^ 


K(i-fo)]cos[v^o-M^„] 

+2 sin  [n^{t- )]  sin  [v^  „  -  „  ] J 


j-2 sin  [n^{t- )] cos [v^ „  - „ ] 

1+  [4 cos  (f  -  f  J]  -  3]  sin  [v^  „  -  , ] J 


-('^c(0-«c)l 


\ sin  [^c  (0  -  44,  (0]  [-3  cos  [n^  (f  -  f„ )]  +  4]  sin  [v,  „  -  „  ] 

[6 sin  [n,  (f  -  )]  -  6n,  (f  -  f„ )]  sin  [v^  „  -  J] 

l+cos[v,„-M^J  J 

fsin[n,(f-f„)] 


+  cosK(0-iWc(0]j 


-(‘ic(0-nc)('ic,o-nc)i 


-sink  (0-44,  (0]l 


-cos[v^„-M^J 


2cos[n,(f-f0]  ^  2 


+  cos[v,(0-M^(0]i 


2cos[n^(f-f0]  2 


4  sin  [«,(?-?„)] 


I  L-3(f-fo) 


0^c,o-44,J 

cos[v,„-M^J 

4nK„-M^J 


O44  =  cos  [v,  (0  - M,  (0]  {cos  [n,  (f  -  )]  cos  [v,  „  - M,  „ ]  +  2 sin  [n,  (f  -  )] sin  [v,  „  - M,  „ ]} 

I  |-2sin  [n,  (f  -  fj]  cos  [v, ,  -  M,,  ] 

1+  [4cos  [n,  (f  -  )]  -  3]  sin  [v,  „  -  M,  J  J 


+  sin 


in[v,(0-M,(0]{ 


+ 


K(0-«c) 


-sin[K,(0-M,(0]{ 


sin[n,(f-f0]  . 

- - ^cos 

V  -M  1 

*^C,0  C,0 

2cos[n, 

)]^  2  " 

■  r  Ji^r  “1 

+ 

—  +  — 

1 

sin  [v,  0- 44,4,  J 

"c 

^c_ 

+  COS 


K(0-44c(0] 


cos 


+ 


2cos[n,(f-f0]  ^ 

-3(f-f0 


4sin[n,(f-fJ] 


[^c,o-44c,o] 

in  |{^c,o  —  44,0  J 
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®45  =cos[vc(f)-M^(f)] 


-cos[?i^(f-fJ]sin[v^„-M^„]  1 
+2 sin  [n^{t-  io )]  cos  [v^  ^  “  ^c.o  ] J 

+  [4 cos  ( f  -  f  J]  -  3]  cos  [v^ ,  -  „  ] 

f  sin[n^(f-fj]  .  p 


-sin[v^(f)-M^(f)]< 


>K(0-^c(0]“ 


b'c.o-^c.o] 


2cos[n^(f-fo)]  2I  r  1 

+ - - ^  +  -  cos[v,,„-M,,] 


«c  f^c] 

4^‘"["°('"")1-3(,-01cos[.,.-M„] 


®  5,  =  -  sin  [^c  (0  -  (0]  [3«c  sin  [n,{t-  f„ )]]  cos  [v^  „  -  „  ] 

+  cos  [v,{t)-M^{t)]  [6n,  cos  [«,  (f  -?„)]-  6nJ  cos  [v^  „  -  „  ] 

-SinK(f)-iWc(0]i  ^  .  r  .  xn  r  if 

4.  X  l+2sin[n,(f-f„)]cos[v,„-M^JJ 

"*"vc,o  ^c)'  r-i  •  r  /  M  •  r  i 

+  COS  V  (t)-M^  (On 

1+  [4 cos  [n,  (f  -  f„ )]  -  3]  cos  [v^  „  -  „  ] 

f  cos  [v,{t)-M^{t)]  [-3  cos  [n^{t-  f„ )]  +  4]  cos  [v,^„  -  ] 


-sin[v,(0-M,(0] 


+  cos[v,(f)-M^(f)] 


-(^c(0-«c)l  .r  .  xnf[6sin[nc(f-fo)]-6«c(f-fo)]cos 

+  sink(0-iWc(0]  r  1 

I  l+sinK„-M^J 


sin[n^(f-fj] 


sinK,o-4^c,o] 


cos[v^(0-M^(f)]< 
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-('^c(0-«c)('^c,o-«c)“ 


2cos[n^(f-f„)]  2 


sinK„-M,,„] 
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cos[v,„-M^J 
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®52  =  sin  [^c  (0  - (0]  [3nc  sin  [«c  )]]  sin  [v,  „  - „ ] 

-  cos  [vc  (0  -  ii^c  (0]  [H  cos  [n^{t- 1, )]  -  6n,  ]  sin  ] 


-('^c,o-«c)l 
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0^3  =-sin[v^(f)-M^(f)] 
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Appendix  F :  MATLAB  Code 


The  numerical  simulations  for  Chapter  V  were  conducted  using  several  original 
MATLAB  programs,  described  below. 

Relative  Satellite  Motion  Plots  and  Comparisons 

The  primary  script  for  calculating  and  comparing  relative  satellite  motion 
trajectories  via  various  models  is  rel_sat_motion5.m.  It  is  the  core  script,  generating  time 
and  true  anomaly  arrays,  calculating  key  constants,  calling  specific  model  functions,  and 
plotting  trajectory  comparisons.  The  function  for  the  two-body  truth  model  is  always 
called.  In  the  input  block,  a  different  switch  may  be  set  to  1  in  order  to  call  the  HCW  or 
Virtual  Chief  models. 

rel_sat_motion5.m  passes  the  following  inputs  to  each  of  the  separate  MATLAB 
functions  described  in  subsequent  sections:  the  chief  satellite  orbit  elements,  the  deputy 
satellite’s  relative  initial  state  in  coordinates  of  the  chiefs  LVLH  frame,  arrays  of  time 
and  chief  true  anomaly,  and  a  few  constant  functions  pre-calculated  by 
rel_sat_motion5  .m. 

Each  of  the  functions  returns  three  arrays,  representing  the  deputy’s  relative 
position  vector  at  each  time  step  in  coordinates  of  the  chiefs  LVLH  frame. 

The  Algorithm. 

This  function  calculates,  plots,  and  compares  relative  satellite  trajectories  from 
multiple  models.  The  algorithm  is  as  follows: 

•  Input.  In  the  input  section,  the  user  must  edit  the  script  in  order  to  select 
which  models  to  call  (Virtual  Chief,  HCW,  or  both),  to  set  the  chiefs  classical  orbit 
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elements,  to  specify  the  deputy's  initial  conditions  in  c-frame  coordinates,  and  to  set  other 
simulation  details  (step  size,  number  of  chief  orbits,  and  epoch  time).  In  the  version  of 
the  code  included  below,  the  input  values  are  designed  to  replicate  Example  n  from 
Shulman  and  Scott  [88]. 

•  Calculations: 

o  Known  constants.  Use  known  two-body  properties  (see  Appendix  D)  to 
compute  the  chiefs  mean  motion,  time  of  perigee  passage,  period,  semilatus 
rectum,  and  angular  momentum. 

o  Newton-Raphson  iteration  to  find  epoch  conditions.  Solve  Kepler's 
equation  (Eq.  (D.4))  to  find  epoch  values  of  chief  eccentric  anomaly,  leading 
to  true  anomaly  and  angular  rate  (via  Eqs.  (D.5)  and  (D.8)). 

o  Discrete  arrays.  Generate  an  evenly  spaced  array  of  chief  true  anomaly 
values,  based  on  the  input  values  for  step  size  and  number  of  chief  orbits.  In 
this  array,  true  anomaly  will  be  allowed  to  take  on  values  greater  than  2n  ,  so 
that  we  can  calculate  and  plot  more  than  one  orbit.  Then  use  two-body 
properties  to  create  arrays  of  chief  eccentric  anomaly,  mean  anomaly,  and 
time,  where  each  element  corresponds  to  an  element  of  the  true-anomaly  array. 
Thus,  for  the  elliptical  chief  problem,  the  time  array  will  not  have  evenly 
spaced  values. 

•  Models.  Organize  variables  into  a  small  number  of  arrays  for  passing  to 
functions.  Call  the  Kepler2.m  function  to  calculate  a  two-body  "truth"  model.  Call  the 
virtual_chief.m  and  HCW.m  functions,  if  selected. 

•  Plots.  Generate  four  figures.  Eirst,  plot  the  c-frame  coordinates  vs.  time 
for  each  selected  model.  Eigure  22  shows  such  a  plot  for  the  example  input  values  below. 
Then,  plot  the  c-frame  coordinates  vs.  chief  true  anomaly,  as  in  Eigure  23. 
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Figure  22.  Coordinates  vs.  Time 


Figure  23.  Coordinates  vs.  Chief  True  Anomaly 

Next,  plot  the  relative  trajectory  projected  into  each  c-frame  coordinate  plane  and  in  three 
dimensions,  as  in  Figure  24. 
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Figure  24.  Relative  Trajectory 

Finally,  plot  the  error  (as  eompared  with  the  two-body  truth  model)  vs.  time  for  the 
Virtual  Chief  and  HCW  models,  as  in  Figure  25.  Also  plot  the  error  magnitude  (the  root 
sum  square  of  the  three  eoordinate  errors)  vs.  time. 
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Figure  25.  Error  vs.  Time 
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The  Code. 


This  is  rel_sat_motion5.m: 


%  Relative  satellite  motion  plotter 

q, 

0 

%  Calculates  and  plots  the  motion  of  a  deputy  satellite  relative  to  a 
%  chief,  according  to  several  different  relative  motion  models.  Each 
model 

%  is  called  by  a  separate  function,  and  the  results  are  plotted  together 
%  versus  time,  versus  chief  true  anomaly,  and  geometrically. 

q, 

o 

%  Input:  Chief  orbit  elements,  deputy  relative  initial  state,  desired 
%  number  of  orbits  and  step  size,  epoch  time 

q, 

0 

%  Units:  SI 


clc;  clear  all;  close  all;  format  long  g 


%  INPUT: 

g, 

0 


%  Desired  models  (set  each  flag  to  1  if  you  wish  to  run  the  model) : 

vc=l;  %  the  Virtual  Chief  method 

hcw=l;  %  Hill ' s-Clohessy-Wiltshire  equations 

%  Chief  orbit  elements: 

a_C=8348000  ;  %m  (chief  semi -major  axis) 

e_C=0.2  ;  %chief  orbit  eccentricity 
i_C=28 . 5*pi/180  ;  %chief  orbit  inclination 
Omega_C=0  ;  %rad  (chief  RAAN) 
omega_C=0  ;  %rad  (chief  argument  of  perigee) 

M_C0=0  ;  %rad  (chief  mean  anomaly  at  epoch) 

%  Deputy  relative  initial  state  (Cartesian  coordinates  in  LVLH  frame 
fixed 

%  to  chief) : 
x0=0  ;  %m 
y0=40000  ;  %m 
z0=40000  ;  %m 
xdot0=0.0111  ;  %m/s 
ydot0=0  ;  %m/s 
zdot0=0  ;  %m/s 

%  Epoch  time: 
t0=0  ;  %s 

%  Desired  number  of  orbits: 

N=3  ; 

%  Desired  step  size,  in  increments  of  chief  true  anomaly: 
ss=pi/360  ;  %rad 
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CALCULATIONS : 


g, 

o 

g, 

0 


%  I .  Known  constants  and  two-body  properties 

mu=3 . 98600441*10''14;  %  (m''3)  /  (s''2)  (gravitational  parameter) 
n_C=sqrt  (mu/ (a_C''3)  )  ;  %rad/s  (chief  mean  motion) 
t_p=tO-M_CO/n_C;  %s  (chief  time  of  perigee  passage) 
P_C=2*pi/n_C;  %s  (chief  period) 

p_C=a_C*  ( l-e_C''2 )  ;  %m  (chief  semilatus  rectum) 
h=sqrt  (mu*p_C)  ;  %(m'^2)/s  (chief  angular  momentum) 


%  II.  Newton-Raphson  iteration  to  find  epoch  conditions 
f lag=0 ; 

E0=M_C0;  %  first  guess 
E1=I0; 

while  flag==0, 

if  abs (El-EO) <10^  (-9) , 
f lag=l ; 
else  E0=E1; 

E1=E0- (E0-e_C*sin (EO) -M_C0) / (l-e_C*cos (EO) ) ; 

end 

end 

E_C0=E1;  %rad  (chief  eccentric  anomaly  at  epoch) 
if  E_C0<0 

E_C0=E_C0+2*pi; 

end 

if  M_C0==0 
E_C0=0; 

end 

s=sin (E_C0) *sqrt (l-e_C^2) / (l-e_C*cos (E_C0) ) ; 
c= (cos (E_C0) -e_C) / (l-e_C*cos (E_C0) ) ; 

nu_C0=atan2 ( s , c) ;  %rad  (chief  true  anomaly  at  epoch) 
if  nu_C0<0 

nu_C0=nu_C0+2  *pi ; 

end 

nudot_C0=n_C*  (l-e_C^2)  ^  (-3/2)  *  (l+e_C*cos  (nu_C0)  )  ''2;  %rad/s  (chief 
angular  rate  at  epoch) 


%  III.  Discrete  arrays 

nu_C=nu_C0 : ss : nu  C0+2*N*pi;  %rad  (chief  true  anomaly  array,  allowed  to 

take  on  values  greater  than  2*pi) 

s2=sin (nu_C) . *sqrt ( l-e_C^2 ) . / ( l+e_C*cos (nu_C) ) ; 

c2= (e_C+cos (nu_C) ) . / ( l+e_C*cos (nu_C) ) ; 

E_C=atan2 ( s2 , c2 ) ;  %rad  (chief  eccentric  anomaly  array) 
i=l;  %atan2  is  defined  from  -pi  to  pi;  I  need  zero  to  2*pi 
while  i<=length (E_C) 
if  E  C (i) <0 
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E_C (i) =2*pi+E_C (i) ; 

end 
i=i+l ; 

end 

M_C=E_C-e_C*sin (E_C) ;  %rad  (chief  mean  anomaly  array) 
t=M_C . /n_C+t_p;  %s  (time  array) 

i2=l;  %I  don't  want  to  reset  time  to  zero  or  negative  after  each  orbit 
while  i2<=N 
i3=2; 

while  i3<=length (t) 
if  t(i3)<t(i3-l) 

t (i3) =t (i3) +P_C; 

end 

i3=i3+l ; 

end 

i2=i2+l; 

end 


%  MODELS: 

g, 

0 


orb_el= [a_C, t_p, i_C, omega_C, e_C, Omega_C] ; 
rel_state= [xO , yO , zO , xdotO , ydotO ,  zdotO ]  ; 
epoch_cond= [tO , nu_C0 , M_C0 , nudot_C0 ] ; 
constants=[h, p_C, mu, n_C, N] ; 

[x, y, z] =Kepler2 (orb_el, rel_state, epoch_cond, t, nu_C, constants) ; 
if  vc==l 

[x_vc, y_vc, z_vc] =virtual_chief (orb_el, rel_state, epoch_cond, t, nu_C, consta 
nts)  ; 

[x_vcc, y_vcc, z_vcc] =vc_check (orb_el, rel_state, epoch_cond, t, nu_C, constant 

s)  ; 

end 

if  hcw==l 

[x_hcw, y_hcw, z_hcw] =HCW (rel_state, epoch_cond, t, constants) ; 

end 


%  PLOTS: 


%  Figure  1  (Coordinates  versus  time) 

subplot ( 3 ,  1 ,  1 )  ;  plot (t, X,  ' k ' ) 

[hand, ob, pi, text] =legend (  '  2BP  '  )  ; 
hold  on; 
if  vc==l 

plot (t, x_vc, ' r-- ' ) 

[hand, ob, pi, text] =legend (hand, text, 'Virtual  Chief ' ) ; 

end 
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if  hcw==l 

plot (t, x_hcw, ' g : ' ) 

[hand, ob, pi , text] =legend (hand, text, ' HCW ' ) ; 

end 

xlabel ( ' t  ( s )  '  ) 
ylabel ( ' x  (m) ' ) 

title (' Relative  coordinates  vs.  time') 
subplot ( 3 , 1 , 2 ) ;  plot (t, y, ' k ' ) 
hold  on; 
if  vc==l 

plot (t, y_vc, ' r-- ' ) 

end 

if  hcw==l 

plot (t, y_hcw, ' g : ' ) 

end 

xlabel ( ' t  ( s )  '  ) 

ylabel ( ' y  (m) ' ) 

subplot ( 3 , 1 , 3 ) ;  plot (t, z, ' k' ) 

hold  on; 

if  vc==l 

plot (t, z_vc, ' r-- ' ) 

end 

if  hcw==l 

plot (t, z_hcw, ' g : ' ) 

end 

xlabel ( ' t  ( s )  '  ) 
ylabel  (  ' z  (m)  ' ) 


%  Figure  2  (Coordinates  versus  chief  true  anomaly) 
figure 

subplot (3,1,1);  plot (nu_C, x, ' k' ) 

[hand2 , ob2 , pl2 , text2 ] =legend ( ' 2BP  '  )  ; 
hold  on; 
if  vc==l 

plot (nu_C, x_vc, ' r-- ' ) 

[hand2 ,  ob2 ,  pl2 ,  text2  ]  =legend  (hand2 ,  text2 ,  'Virtual  Chief) 

end 

if  hcw==l 

plot (nu_C, x_hcw, ' g : ' ) 

[hand2  ,  ob2  ,  pl2 , text2 ] =legend (hand2 , text2 ,  ' HCW ' ) ; 

end 

title (' Relative  coordinates  vs.  chief  true  anomaly') 
xlabel ( ' nu_C  (rad) ' ) 
ylabel ( ' x  (m) ' ) 

subplot (3,1,2);  plot (nu_C, y, ' k' ) 
hold  on; 
if  vc==l 

plot (nu_C, y_vc, ' r-- ' ) 

end 

if  hcw==l 

plot (nu_C, y_hcw, ' g : ' ) 

end 

xlabel ( ' nu_C  (rad) ' ) 
ylabel ( ' y  (m) ' ) 

subplot (3,1,3) ;  plot (nu_C, z, ' k' ) 
hold  on; 
if  vc==l 
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plot (nu_C, z_vc, ' r-- ' ) 

end 

if  hcw==l 

plot (nu_C, z_hcw, ' g ; ' ) 

end 

xlabel ( ' nu_C  (rad) ' ) 
ylabel ( ' z  (m) ' ) 


%  Figure  3  (Geometric  plots) 
figure 

subplot (2 , 2 , 1 ) ;  plot (x,  y,  ' k ' ) 
axis  equal 
xlabel ( ' X  (m) ' ) 
ylabel ( ' y  (m) ' ) 

[hand3 , ob3 , pl3 , text3 ] =legend ( ' 2BP ' ) ; 
hold  on; 
if  vc==l 

plot (x_vc, y_vc, ' r-- ' ) 

[hand3 , ob3 , pl3 , text3 ] =legend (hand3 , text3 , 'Virtual  Chief) 

end 

if  hcw==l 

plot (x  hew, y_hcw, ' g ; ' ) 

[hand3 , ob3 , pl3 , text3 ] =legend (hand3 , text3 , ' HCW ' ) ; 

end 

plot (x (1) , y  (1) ,  '  ko  '  ) 
plot(x(15) ,y(15) ,  'kx') 
if  vc==l 

plot (x_vc ( 15 ) , y_vc ( 15 ) , ' rx ' ) 

end 

if  hcw==l 

plot (x_hcw ( 15 ) , y_hcw ( 15 )  ,  ' gx '  ) 

end 

title (' Relative  trajectory') 

subplot (2 , 2 , 2 ) ;  plot (x,  z ,  ' k ' ) 

axis  equal 

xlabel ( ' X  (m) ' ) 

ylabel ( ' z  (m) ' ) 

hold  on; 

if  vc==l 

plot (x_vc, z_vc, ' r-- ' ) 

plot (x_vc ( 15 ) , z_vc ( 15 ) , ' rx ' ) 

end 

if  hcw==l 

plot (x  hew, z_hcw, ' g : ' ) 

plot (x_hcw ( 15 ) , z_hcw ( 15 ) , ' gx ' ) 

end 

plot (x (1) , z  (1)  ,  '  ko  '  ) 
plot(x(15) ,z(15) ,  'kx') 
subplot (2, 2, 3) ;  plot (y, z, ' k' ) 
axis  equal 
xlabel ( ' y  (m) ' ) 
ylabel ( ' z  (m) ' ) 
hold  on; 
if  vc==l 

plot (y_vc, z_vc, ' r-- ' ) 

plot (y_vc ( 15 ) , z_vc ( 15 ) , ' rx ' ) 

end 
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if  hcw==l 

plot (y_hcw, z_hcw,  ' g :  ' ) 

plot (y_hcw ( 15 ) , z_hcw ( 15 )  ,  ' gx  '  ) 

end 

plot (y (1) , z  (1)  ,  '  ko  '  ) 

plot(y(15) ,z(15) ,  'kx') 

subplot  (2,2,4);  plot3(x,y,z,'k') 

xlabel ( ' X  (m) ' ) 

ylabel ( ' y  (m) ' ) 

zlabel ( ' z  (m) ' ) 

hold  on; 

if  vc==l 

plots (x_vc, y_vc, z_vc, ' r-- ' ) 

plots (x_vc ( 15 ) , y_vc ( 15 ) , z_vc ( 15 )  ,  ' rx ' ) 

end 

if  hcw==l 

plots (x_hcw, y_hcw, z_hcw, ' g ; ' ) 

plots (x_hcw ( 15 ) , y_hcw ( 15 ) , z_hcw ( 15 )  ,  ' gx ' ) 

end 

plots  (x  (1)  ,  y  (1)  ,  z  (1)  ,  '  ko  '  ) 
plots (x(15) ,y(15) , z (15)  ,  'kx'  ) 
axis  equal 
grid  on 


%  Figure  4  (Error) 
figure 

subplot  (4,1,1) 
plot ( 0 , 0 ,  ' w  ' ) 
hold  on 

[hand4 , ob4 , pl4 , text4 ] =legend ( '  ' ) ; 

if  vc==l 

plot (t, x_vc-x, ' r-- ' ) 

[hand4 ,  ob4 ,  pl4 ,  text4  ]  =legend  (hand4 ,  text4 ,  'Virtual  Chief) 

end 

if  hcw==l 

plot (t, x_hcw-x, ' g : ' ) 

[hand4 , ob4 , pl4 , text4 ] =legend (hand4 , text4 , ' HCW ' ) ; 

end 

legend  boxoff 
xlabel ( ' t  ( s )  '  ) 
ylabel ('x  error  (m) ') 

title ('Error  as  compared  with  two-body  solution') 
subplot ( 4 , 1,2) 
hold  on 
if  vc==l 

plot (t, y_vc-y, ' r-- ' ) 

end 

if  hcw==l 

plot (t, y_hcw-y, ' g : ' ) 

end 

xlabel ( ' t  ( s )  '  ) 
ylabel ('y  error  (m) ') 
subplot ( 4 , 1 , S ) 
hold  on 
if  vc==l 

plot (t, z_vc-z, ' r-- ' ) 

end 
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if  hcw==l 

plot ( t , z_hcw-z , ' g : ' ) 

end 

xlabel ( ' t  ( s )  '  ) 
ylabel  (  ' z  error  (m)  ' ) 
subplot (4,1,4) 
hold  on 
if  vc==l 

err_vc=sqrt  (  (x_vc-x)  .  ''2+  (y_vc-y)  .  ''2+  (z_vc-z)  .  ''2)  ; 
plot (t, err_vc, ' r-- ' ) 
rms_vc=sqrt  (mean  (err_vc  .''2)) 

end 

if  hcw==l 

err_hcw=sqrt  (  (x_hcw-x)  .  ''2+  (y_hcw-y)  .  ''2+  (z_hcw-z)  .  ''2)  ; 
plot (t, err_hcw, ' g : ' ) 
rms_hcw=sqrt  (mean  (err_hcw .''2)) 

end 

xlabel ( ' t  ( s )  '  ) 

ylabel (' error  magnitude  (m) ') 


The  Two-body  Keplerian  Truth  Model 

The  Algorithm. 

This  function  calculates  the  relative  satellite  trajeetory  using  the  exaet  solution  to 
the  nonlinear  two-body  problem.  The  algorithm  is  as  follows: 

•  1.  Calculate  the  rotation  matrix  from  the  chief  s  LVLH  frame  to  the  i- 

frame  (Earth-eentered  inertial),  evaluated  at  epoeh.  See  below  for  derivation  of  the 
matrix. 


•  n.  Using  the  matrix  just  calculated,  compute  the  chiefs  state  vectors  in  i- 
frame  coordinates: 


^c/ 


=  R. 


^C/ 


dt  ^/i 


=  R. 


'  d  ^ 


•  111.  Using  the  same  matrix,  compute  the  deputy’s  relative  state  vectors  in 

i-frame  coordinates: 


145 


V  =^c  V 

/c  ji  L  /c 


d/% 


rv.  Add  these  veetors  to  find  the  deputy’s  inertial  state: 


%L“  Vc 


d  ^  '  d  ^ 


—  f)/  =  — f’/  +  — A)/ 
dt  //  dt  ^  dt  /c 


•  V.  Use  known  two-body  properties  (see  Appendix  D)  to  calculate  the 
deputy’s  orbit  elements,  ,  co^ ,  and  q  ,  as  well  as  related  values  £^  0 . 


tpj,,  and  . 


•  VI.  Using  the  time  array  passed  from  rel_sat_motion.m,  generate  an  array 
of  values  for  the  deputy’s  mean  anomaly: 


Md  (t)  njj  ^ 


•  Vn.  Using  Newton-Raphson  iteration,  compute  the  deputy’s  eccentric  and 
true  anomalies  at  each  step: 


^  oif) 

{A  sin  E^{t) 

sm  v„  U  = 

1-e^  cosEjj  (t) 


:os£^(t)-g^ 
-  e cos 


Vni.  At  each  time  step,  compute  the  chief  and  deputy  orbit  radius: 


Pc 

+  e^  cosUp  (0 


,M  = _ Pp _ 

1-i-g^  cosv^  (t) 
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•  IX.  At  each  time  step,  calculate  a  rotation  matrix  from  LVLH  to  i-frame 
coordinates,  evaluated  at  the  current  time,  for  both  the  chief  and  the  deputy.  This  matrix 
follows  the  same  formula  as  that  in  Step  I. 

•  X.  At  each  time  step,  rotate  the  chief  and  deputy  position  vectors  into  i- 
frame  coordinates,  as  was  done  for  the  chief  initial  position  in  Step  n. 

•  XL  At  each  time  step,  difference  the  chief  and  deputy  positions  to  create  a 
relative  position  vector  in  i-frame  coordinates: 


•  XII.  At  each  time  step,  use  the  transpose  of  the  chiefs  rotation  matrix 
from  Step  IX  to  rotate  the  relative  position  vector  into  coordinates  of  the  chiefs  LVLH 
frame: 


x{t) 

y(') 


The  Rotation  Matrices. 

The  rotation-matrix  formula  used  in  Steps  I  and  IX  of  the  algorithm  will  now  be 
derived. 

An  arbitrary  vector  r  defined  in  i-frame  (ECI)  coordinates  can  be  expressed  in 
coordinates  of  a  particular  satellite’s  nodal-equatorial  (q)  reference  frame  as  follows: 


[7] 


where  R'^‘’  is  simply  a  three -rotation  of  magnitude  Q ,  the  satellite’s  right  ascension  of 
the  ascending  node: 
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R‘^'^  =  (Q) 


cos  Q  sin  Q  0 
-sinQ  cosQ  0 
0  0  1 


(F.l) 


This  rotation  is  illustrated  in  Figure  26,  which  depicts  the  i-frame  unit  basis  vector 
ij  rotated  to  the  q-frame. 


1 

0 

0 


(sinQ)^2 


Figure  26.  Rotation  from  i  to  q 

Likewise,  a  vector  defined  in  q-frame  coordinates  can  be  expressed  in  the 
satellite’s  nodal  (n)  frame  via  a  1 -rotation  of  magnitude  i ,  the  satellite’s  inclination; 
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j^q^n 


(F.2) 


=  ^i(0 


1  0  0 
0  cos  i  sin  i 
0  -sinz  cosz 


This  rotation  is  illustrated  in  Figure  27,  which  depicts  the  q-frame  unit  basis 
vector  rotated  to  the  n-frame. 


n 

yj 


y' 


0 

1 

0 


=  (cos  z )  /22  “  (sin  /) 


Figure  27.  Rotation  from  q  to  n 

It  is  worth  noting  that  for  satellites  in  equatorial  orbits,  Q  and  z  should  each  be 
set  to  zero.  The  rotations  in  Eqs.  (F.l)  and  (F.2)  will  then  produce  the  identity  matrix. 

Finally,  a  vector  defined  in  n-frame  coordinates  can  be  expressed  in  the  satellite’s 
orbit  (o)  frame  via  a  3-rotation  of  magnitude  u  ,  the  satellite’s  argument  of  latitude.  If  the 
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satellite  is  in  a  non-circular  orbit,  u  may  be  replaced  by  its  equivalent,  co  +  v ,  ox 
argument  of  perigee  plus  true  anomaly: 


=R^{o)  +  v) 


COs(<i)-|- v) 
-sin(n;-l- v) 
0 


sin  (ft) -I- v) 
cos(ft)-l- v) 
0 


0 

0 

1 


(F.3) 


Thus,  we  can  use  o-frame  coordinates  to  express  a  vector  defined  in  i-frame 
coordinates  as  follows: 


=  7?'^  [r],  =  [7]^ 


Multiplying  the  matrices  from  Eqs.  (F.l),  (F.2),  and  (F.3)  yields 


R 


i—>o 


cQc(ft)-i- v)-sQcis(ft)-i- v) 
-cQs(ft)-i- v)-sQc/c(ft)-i- v) 
sQsi 


s  Qc  (ft)-i- v) -I- cQc  /  s  (ft)-i- v) 
-  s  Q  s  ( ft) -I- V ) -I- c  Q  c  /  c  ( ft) -I- V ) 
-cQs/ 


sis(ft)-i- v) 
sic(ft)-i- v) 
ci 


where  cosine  and  sine  are  signified  by  c  and  s,  respectively.  This  is  the  rotation  needed 
for  the  chief  satellite  in  Step  XII  of  the  algorithm. 

However,  we  need  the  inverse  transformation,  R“^‘ .  Since  the  inverse  of  a 
rotation  matrix  is  its  transpose. 


R 


o^i 


cQc(ft)-i-v)-sQc/s(ft)-i- v) 
sQc(ft)-i- v)  +  cQc/s(ft)-i-v) 
sis[co  +  v) 


-cQs(ft)-i-v)-sQc/c(ft)-i- v) 
-  s  Q  s  ( ft) -I- V ) -I- c  Q  c  /  c  ( ft) -I- V ) 
sic[co  +  v) 


sQsi 

-cQsi 

ci 


(F.4) 
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Eq.  (F.4)  gives  the  formula  needed  in  Steps  I  and  IX  of  the  algorithm. 


The  Code. 

This  is  Kepler2.m,  the  MATLAB  function  which  calculates  the  relative  satellite 
trajectory  using  fully  nonlinear  two-body  Keplerian  motion: 

%  Two-body  Exact  Relative  Orbit  Plotter  (rev  2) 

"5 

%  Provides  the  exact  relative  trajectory  of  a  deputy  satellite  in  the 
%  local-vertical/local-horizontal  (LVLH)  reference  frame  of  a  chief 
%  satellite,  assuming  simple  two-body  Keplerian  motion. 

%  This  function  should  work  for  all  chief  orbits  (circular  or  elliptical, 

%  equatorial  or  inclined) . 

'6 

%  Required  inputs  are  the  chief's  orbit  elements,  the  deputy's  relative 
%  initial  state  in  LVLH  coordinates,  and  matched  arrays  of  time  and 
chief 

%  true  anomaly. 

function  [x, y, z] =Kepler2 (orb_el, rel_state, epoch_cond, t, nu_C, constants) 

%  Assign  variables 

a_C=orb_el (1 ) ; 
i_C=orb_el (3) ; 
omega_C=orb_el (4) ; 
e_C=orb_el (5) ; 

Omega_C=orb_el ( 6 ) ; 
xO=rel_state (1 ) ; 
yO=rel_state (2 ) ; 
zO=rel_state (3) ; 
xdotO=rel_state ( 4 ) ; 
ydotO=rel_state (5) ; 
zdotO=rel_state ( 6) ; 
tO=epoch_cond ( 1 ) ; 
nu_CO=epoch_cond (2 ) ; 
nudot_CO=epoch_cond ( 4 ) ; 
p_C=constants (2) ; 
mu=constants (3) ; 


%  I.  Rotation  matrix  from  chief  LVLH  frame  (o)  to  ECI  frame  (i) ,  evaluated 
%  at  epoch.  Uses  R_C=R(o  to  i)=R' (i  to  o) , 

%  where  R(i  to  o ) =R3 (omega+nu) R1 ( i ) R3 (Omega) . 

R_C0 (1,1) =cos (Omega_C) *cos (omega_C+nu_CO ) - 
sin (Omega_C) *cos (i_C) *sin ( omega_C+nu_CO ) ; 

R_C0 (1,2) =-cos (Omega_C) *sin ( omega_C+nu_CO ) - 
sin (Omega_C) *cos (i_C) *cos ( omega_C+nu_CO ) ; 

R_C0 (1 , 3 ) =sin (Omega_C) *sin (i_C) ; 

R_C0  (2 , 1 ) =sin (Omega_C) *cos (omega_C+nu_CO ) +cos (Omega_C) *cos (i_C) *sin ( omega_C+nu_C 

0)  ; 

R_C0 (2,2)=- 

sin (Omega_C) *sin (omega_C+nu_C0 ) +cos (Omega_C) *cos (i_C) *cos ( omega_C+nu_C0 ) ; 

R_C0 (2 , 3 ) =-cos (Omega_C) *sin (i_C) ; 

R_C0 (3, 1 ) =sin (i_C) *sin (omega_C+nu_C0 ) ; 

R_C0 (3 , 2 ) =sin (i_C) *cos (omega_C+nu_C0 ) ; 

R  CO (3, 3) =cos  (i  C)  ; 
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%  II.  Calculate  chief  epoch  state  in  Earth-centered  inertial  (ECI) 

%  coordinates 

r_C0=p_C/ ( l+e_C*cos (nu_C0 ) ) ;  %m  (chief  orbit  radius  at  epoch) 
pos_C0_o= [r_C0 ; 0 ; 0 ] ;  %  chief  position  vector  at  epoch,  LVLH  coordinates 
rdot_CO=sqrt (inu/p_C) *e_C*sin (nu_C0 ) ;  %m/s  (chief  radial  rate  at  epoch) 
vel_C0_o= [rdot_C0 ; r_C0*nudot_C0; 0 ] ;  %  chief  velocity  vector  at  epoch,  LVLH 
coordinates 

pos_C0_i=R_C0*pos_C0_o;  %  chief  position  vector  at  epoch,  ECI  coordinates 
vel_C0_i=R_C0*vel_C0_o;  %  chief  velocity  vector  at  epoch,  ECI  coordinates 


%  III.  Rotate  deputy  relative  state  to  ECI  coordinates 

pos_relO_i=R_CO* [xO ; yO ; zO ] ;  %  relative  position  vector  at  epoch,  ECI  coordinates 
vel_relO_i=R_CO* ( [xdotO; ydotO; zdotO] +cross ( [ 0 ; 0; nudot_C0 ] , [xO ; yO ; zO ] ) ) ;  % 
relative  velocity  vector  at  epoch,  ECI  coordinates 


%  IV.  Add  vectors  to  find  deputy  inertial  state 

pos_D0_i=pos_C0_i+pos_rel0_i ;  %  deputy  position  vector  at  epoch,  ECI  coordinates 
vel_D0_i=vel_C0_i+vel_rel0_i ;  %  deputy  velocity  vector  at  epoch,  ECI  coordinates 


%  V.  Calculate  deputy  orbit  elements 

r_D0=norm (pos_D0_i ) ;  %m  (deputy  radius  at  epoch) 
v_D0=norm (vel_D0_i ) ;  %m/s  (deputy  speed  at  epoch) 
epsilon_D=v_D0^2/2-mu/r_D0 ;  %m^2/s^2  (deputy  specific  energy) 
a_D=-mu/ (2*epsilon_D)  %m  (deputy  semi-major  axis) 

H_D_vec=cross (pos_D0_i , vel_D0_i ) ;  %deputy  angular  momentum  vector 
H_D=norm (H_D_vec ) ;  %m^2/s  (deputy  angular  momentum) 

e_D_vec=cross ( vel_D0_i , H_D_vec ) /mu-pos_D0_i/r_D0 ;  %  deputy  eccentricity  vector 
e_D=norm (e_D_vec )  %deputy  eccentricity 

i_D=acos(([0  0  1 ] *H_D_vec) /H_D)  %rad  (deputy  inclination) 

nhat=cross ([ 0 ; 0; 1 ], H_D_vec) /norm (cross ([ 0 ; 0; 1 ], H_D_vec) ) ;  %deputy  unit  nodal 
vector 

Omega_D=atan2 (nhat (2 ) , nhat (1 ) )  %rad  (deputy  RAAN) 

omega_D=acos ( (e_D_vec/e_D) ' *nhat)  %rad  (deputy  argument  of  perigee) 
if  [0  0  l]*e_D_vec<0 

omega_D=2*pi-omega_D 

end 

nu_D0=acos ( (e_D_vec/e_D) ' * (pos_D0_i/r_D0 ) )  %rad  (deputy  true  anomaly  at  epoch) 
if  pos_D0_i ' *vel_D0_i<0 

nu_D0=2*pi-nu_D0  %  correcting  for  cases  when  the  deputy  is  past  apogee 

end 

E_D0=2*atan2 (sqrt ( (l-e_D) / (l+e_D) ) *sin (nu_D0/2) , cos (nu_D0/2) ) ;  %rad  (deputy 
eccentric  anomaly  at  epoch) 

n_D=sqrt (mu/a_D^3) ;  %rad/s  (deputy  mean  motion) 

t_pD=tO- (E_D0-e_D*sin (E_D0 ) ) /n_D;  %s  (deputy  time  of  perigee  passage) 
p_D=a_D* (l-e_D*2 ) ;  %m  (deputy  semilatus  rectum) 


%  VI .  Create  deputy  mean  anomaly  array 

M_D=n_D. * (t-t_pD) ;  %  rad  (deputy  mean  anomaly  array) 
for  i=l : length (M_D)  %I  want  to  reset  M_D  to  zero  after  each  orbit 
f=fix (M_D (i) / (2*pi) ) ; 

M_D (i) =M_D (i) -2*f*pi; 

end 

%  VII.  Create  deputy  true  anomaly  array 
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for  12=1 : length (t) 
flag=0 ; 

E0=M_D(i2);  %  first  guess 
El=8; 

while  flag==0, 

if  abs  (El-E0)<10^ (-9) , 
flag=l ; 
else  E0=E1; 

E1=E0- (E0-e_D*sin (EO) -M_D (12) ) / (l-e_D*cos (EO) ) ; 

end 

end 

E_D(i2)=El;  %rad  (deputy  eccentric  anomaly) 
s2=sin (E_D (12 )  ) *  sqrt ( l-e_D^2 ) / ( l-e_D*cos (E_D (12))); 
c2= (cos (E_D (12) ) -e_D) / (l-e_D*cos (E^D (i2) ) )7 
nu_D (i2) =atan2 (s2, c2) ;  %rad  (deputy  true  anomaly) 

end 


%  Propagate  chief  and  deputy  inertial  states,  then  rotate  to  LVLH 
for  i3=l : length (t) 

%  VIII.  Chief  and  deputy  trajectory  and  velocity  equations 
r_C=p_C/ (lte_C*cos (nu_C (13) ) ) ;  %m  (chief  orbit  radius) 
r_D=p_D/ (l+e_D*cos (nu_D (13) ) ) ;  %m  (deputy  orbit  radius) 
v_C=sqrt (2*mu/r_C-mu/a_C) ;  %m/s  (chief  velocity) 
v_D=sqrt (2*mu/r_D-mu/a_D) ;  %m/s  (deputy  velocity) 

%  IX.  Chief  and  deputy  rotation  matrices 

%  rotation  matrix  from  chief  LVLH  frame  to  ECI  frame,  using  R(i  to 
o) =R3 (omega+nu) R1 (i) R3 (Omega)  and  R(o  to  i)=R'(i  to  o) 

R_C (1,1) =cos (Omega_C) *cos ( omega_C+nu_C (13) ) - 
sin (Omega_C) *cos (i_C) *sin (omega_C+nu_C (13) ) ; 

R_C (1,2) =-cos (Omega_C) *sin (omega_Ctnu_C (13) ) - 
sin (Omega_C) *cos (i_C) *cos (omega_C+nu_C (13) ) ; 

R_C (1,3) =sin (Omega_C) *sin (i_C) ; 

R_C (2,1) =sin (Omega_C) *cos (omega_C+nu_C (13) ) tcos (Omega_C) *cos (i_C) *sin (omega_C+nu 
_C(i3) ) ; 

R_C (2,2)=- 

sin (Omega_C) *sin (omega_C+nu_C (i3) ) tcos (Omega_C) *cos (i_C) *cos (omega_C+nu_C (13) ) ; 
R_C (2, 3) =-cos (Omega_C) *sin  (i_C) ; 

R_C (3,1) =sin (i_C) *sin (omega_C+nu_C (13) ) ; 

R_C (3, 2) =sin (i_C) *cos (omega_C+nu_C (13) ) ; 

R_C  (3,  3)  =cos  (i_C)  ; 

%  rotation  matrix  from  deputy  orbit  frame  to  ECI  frame 
R_D (1,1) =cos (Omega_D) *cos ( omega_D+nu_D (13) ) - 
sin (Omega_D) *cos (i_D) *sin (omega_D+nu_D (13) ) ; 

R_D ( 1 , 2 ) =-cos (Omega_D) *sin (omega_D+nu_D ( 13 ) ) - 
sin (Omega_D) *cos (i_D) *cos (omega_D+nu_D (13) ) ; 

R_D ( 1 , 3 ) =sin (Omega_D) *  sin ( i_D) ; 

R_D ( 2 , 1 ) =sin (Omega_D) *cos (omega_D+nu_D (13 ) ) tcos (Omega_D) *cos (i_D) *  sin ( omega_Dtnu 
_D(i3)); 

R_D (2,2)=- 

sin (Omega_D) *sin (omega_Dtnu_D ( 13 ) ) tcos (Omega_D) *cos ( i_D) *cos (omega_Dtnu_D ( 13 ) )  ; 
R_D (2, 3) =-cos (Omega_D) *sin  (i_D) ; 

R_D ( 3, 1 ) =sin (i_D) *  sin (omega_Dtnu_D (13 ) ) ; 

R_D (3, 2) =sin (i_D) *cos (omega_Dtnu_D (13) ) ; 

R_D (3,3) =cos (i_D) ; 

%  X.  Chief  and  deputy  conditions  rotated  to  ECI  frame 
r_C_vec=R_C* [r_C ; 0 ; 0 ] ;  %  chief  ECI  position  vector 
r_D_vec=R_D* [r_D; 0 ; 0 ] ;  %  deputy  ECI  position  vector 

%  XI .  ECI  relative  position  vector 
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r_EC  I_ve  c=r_D_ve  c-  r_C_ve  c  ; 

%  XII.  Relative  coordinates  rotated  back  to  chief  LVLH  frame 

r_LVLH_vec=R_C ' *r_ECI_vec;  %  relative  position  vector  in  chief  orbit  frame 
coordinates 

X ( i3 ) =r_LVLH_vec ( I ) ;  %m 
y  ( i3 ) =r_LVLH_vec (2 )  ;  %m 
z  ( i3 ) =r_LVLH_vec (3 )  ;  %m 

end 

Examples  Compared  against  Satellite  Toolkit. 

In  order  to  demonstrate  the  accuracy  of  the  two-body  relative  satellite  motion 
predicted  by  Kepler2.m,  I  will  compare  its  output  against  output  from  the  Two  Body 


propagator  in  Satellite  Toolkit  (STK).  I  will  perform  this  comparison  for  three  examples. 


The  first  example  is  a  typical  relative  satellite  motion  case,  with  small  distances 


and  small  but  non-negligible  eccentricities.  For  this  example,  the  orbit  elements  of  the 


chief  satellite  and  the  deputy  satellite’s  relative  state  in  coordinates  of  the  chief  LVLH 


frame  are  given  in  Table  7,  along  with  the  resulting  deputy  orbit  elements. 


Table  7.  Input  for  the  First  Example 


Chief  orbit 

Relative  initial  state 

Deputy  orbit 

Of.  =  26778137  m 

Xo= -2357.02260395516  m 

=  26778090.7194924  m 

-  0.01 

yo=  5714.04520791032  m 

=0.0100867011056697 

ic  =  28.5° 

o 

II 

o 

=0.49756671315498  rad 

Q^=0 

Xo  =0.35626933756075  — 

s 

=6.67858183316407e-008  rad 

II 

o 

yo  =0.686069106910399  — 

s 

co^  =  6.27424251721299  rad 

^c,o=0 

Zo  =0.576312899024239  — 

s 

0  =  0.0091582905573582  rad 
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Results  for  the  first  example  are  compared  over  four  chief  orbits  and  shown  in 
Figure  28.  For  this  simulation,  the  epoch  time  was  set  to  zero,  and  the  chief  true  anomaly 
step  size  was  set  to  half  a  degree. 


The  second  example  illustrates  a  case  with  large  separation  distances  and  a  chief 
in  a  circular,  equatorial  orbit.  For  this  example.  Table  8  gives  the  orbit  elements  and 
initial  relative  coordinates. 

Results  for  the  second  example  are  compared  over  four  chief  orbits  and  shown  in 
Figure  29.  Epoch  time  is  zero,  and  chief  true  anomaly  step  size  is  half  a  degree. 
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Table  8.  Input  for  the  Second  Example 


Chief  orbit 

Relative  initial  state 

Deputy  orbit 

=  6778137  m 

Xq  =  2000  m 

=679031 1.93490504  m 

Cc  =0 

=100000  m 

=0.00139062906315371 

o 

II 

Zo  =  2000 

ij,  =  0.000304358514095287  rad 

Xq  =0.35626933756075  — 

s 

=  -1.30671940634345  rad 

II 

o 

yo  =0.686069106910399  — 

s 

0^  =1.28706041049096  rad 

Mc.o=0 

Zo  =0.576312899024239  — 

s 

v^o  =0-0344069021226111  rad 
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The  third  example  illustrates  satellites  in  highly  eccentric  orbits.  Table  9  gives  the 


orbit  elements  and  initial  relative  coordinates. 


Table  9.  Input  for  Third  Example 


Chief  orbit 

Relative  initial  state 

Deputy  orbit 

=  30778137  m 

Xg  =  0  m 

=30777601. 1837545  m 

-  0.75 

Jq  =150  m 

=0.749999394965603 

=28.5° 

Zq  =  2000  m 

=  0.497453372017292  rad 

o 

II 

=  0.000390591605232089  rad 

II 

o 

o 

II 

o 

co^  =  6.28283388440474  rad 

Mf.  n  =  —  rad 

C.O  g 

m 

Zo  =1  — 
s 

o  =1.97382718065585  rad 

Results  for  the  third  example  are  compared  over  three  chief  orbits  and  shown  in 
Figure  30.  Again,  epoch  time  is  zero,  and  chief  true  anomaly  step  size  is  half  a  degree. 

As  demonstrated  by  these  three  widely  varying  examples,  the  relative  trajectories 
calculated  by  Kepler2.m  provide  an  accurate  representation  of  the  nonlinear  two-body 
solution. 


157 


The  Virtual  Chief  Method 


The  Algorithm. 

This  function  calculates  the  relative  satellite  trajectory  using  the  Virtual  Chief 
model.  The  algorithm  is  as  follows: 

•  1.  Calculate  the  intermediate  parameters  A,  B  ,  C  ,  and  D  .  Then  print  a 

boundedness  check:  if  motion  is  bounded,  then  the  index  6 A  +  3C  should  equal  zero. 
Next,  calculate  and  print  the  virtual  chief  parameters  A, ,  (j)^,  A^,  (f)^,  ,  and  . 

•  n.  Compute  the  solution,  the  deputy's  c-frame  position  coordinates,  for 
each  time  step. 
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The  Code. 


This  is  virtual_chief.m: 


%  Virtual  chief  model  relative  satellite  motion  plotter 

q, 

0 

%  Assumptions:  Both  satellites  remain  close  to  virtual  chief  in  a 
circular 

%  orbit  (i.e.,  small  separation  and  small  eccentricity) . 

q, 

0 

%  Model  initial  time  is  t=tO . 


%  Reference  frame:  The  c-frame,  a  Cartesian  LVLH  frame  attached  to 
actual 

%  chief.  Realization  is  the  six  new  parameters. 


function 

[x, y, z] =virtual_chief (orb_el, rel_state, epoch_cond, t, nu_C,  constants) 

%  Assign  Variables 

t_p=orb_el (2) ; 
e_C=orb_el (5) ; 
xO=rel_state (1) ; 
yO=rel_state  (2) ; 
zO=rel_state  (3) ; 
xdotO=rel_state  (4) ; 
ydotO=rel_state  (5) ; 
zdotO=rel_state (6) ; 
tO=epoch_cond ( 1 ) ; 
nu_CO=epoch_cond (2 ) ; 

M_CO=epoch_cond ( 3 ) ; 
nudot_CO=epoch_cond ( 4 ) ; 
n  C=constants  ( 4 ) ; 


%  I.  Parameters: 

A=cos (nu_C0-M_C0 ) *xO-sin (nu_C0-M_C0 ) *y0 ; 

B=sin (nu  CO-M  CO ) *xOtcos (nu  C0-M_C0 ) *y0 ; 

C= (l/n_C) * (sin (nu_C0-M_C0) * (xdotO- (nudot_C0-n_C) *y0) +cos (nu_C0- 
M_C0) * (ydotO+ (nudot_C0-n_C) *x0) ) ; 

D= (l/n_C) * (cos (nu_C0-M_C0) * (xdotO- (nudot_C0-n_C) *y0) -sin (nu_C0- 
M_C0 ) * (ydotO+ (nudot_C0-n_C) *x0 ) ) ; 

6*A+3*C  %boundedness  check 

Al=sqrt  (  (9/4)  *D'^2+ (81/4)  *A^2+27*A*Ct9*C'^2 ) 
phil=n_C*t0tatan2 (3*A+2*C, D) 

A2=sqrt  (B'^2-4*B*D+4*D'^2  +  16*A''2  +  16*A*Ct4*C'^2) 
phi2=atan2 (4*A+2*C, B-2*D) 
z_max=sqrt  (z0'^2  +  (zdotO/n_C)  ''2) 
psi0=atan2 (zO, zdotO/n_C) 
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Solution : 


%  II. 


x=AI*sin (nu_C+n_C*t_p-phil ) +A2*sin (nu_C- 

n_C*t+n_C*t_p+phi2 ) + (AI/3) *sin (nu_C-2*n_C*t+n_C*t_p+phil ) - 
(3/2) *A2*sin (phi2) *n_C* (t-tO) . *sin (nu_C-n_C*t+n_C*t_p) ; 
y=AI*cos (nu_C+n_C*t_p-phil ) +A2*cos (nu_C- 

n_C*t+n_C*t_p+phi2 ) + (AI/3) *cos (nu_C-2*n_C*t+n_C*t_p+phil ) - 
(3/2) *A2*sin (phi2) *n_C* (t-tO) . *cos (nu_C-n_C*t+n_C*t_p) ; 
z=z  max*sin (n_C* (t-tO ) +psiO )  ; 


Hill-Clohessv-Wiltshire  equations 

The  Algorithm. 

This  function  calculates  the  relative  satellite  trajectory  using  the  HCW  model. 
The  algorithm  is  as  follows: 

•  1.  Calculate  the  submatrices  of  the  state  transition  matrix,  defined  in  Eqs. 

(3.4)  and  (3.6). 


•  n.  Using  the  submatrices,  compute  the  solution  for  each  element  in  the 
discrete  time  array.  The  result  will  be  the  deputy's  relative  state  vectors  in  c-frame 
coordinates  for  each  time  step. 


•  ni.  Calculate  and  print  the  values  of  the  Relative  Orbit  Elements 
described  in  Chapter  III. 

The  Code. 


This  is  HCW.m: 


%  Clohessy-Wiltshire  Relative  Satellite  Motion  Plotter 

g, 

0 

%  Based  on  Equations  3.34  -  3.39  from  Wiesel,  _Spaceflight  Dynamics_, 

2nd 
%  Ed. 

g, 

0 

%  Assumptions:  Chief  satellite  in  circular  orbit;  distance  from  chief  to 
%  deputy  satellite  is  small  compared  to  distance  to  center  of  the  Earth. 

g, 

0 

%  Reference  frame:  converted  to  Cartesian  coordinates  fixed  to  the  chief 
(radial, 

%  along-track,  and  normal  to  the  orbit  plane) . 

g, 

0 

%  Units:  converted  to  SI.  After  conversions,  these  results  should  match 
%  Algorithm  47  from  Vallado,  _Fundamentals  of  Astrodynamics  and 
%  Applications_,  2nd  Ed. 
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function  [x, y, z] =HCW (rel_state, epoch_cond, t, constants) 

%  Assign  variables 

xO=rel_state (1) ; 
yO=rel_state  (2) ; 
zO=rel_state  (3) ; 
xdotO=rel_state  (4) ; 
ydotO=rel_state  (5) ; 
zdotO=rel_state (6) ; 
tO=epoch_cond ( 1 ) ; 

n=constants ( 4 ) ;  %  chief  mean  motion 
%Deputy  initial  relative  position: 
delta_rO= [xO ; yO ; zO ] ;  %m 
%Deputy  initial  relative  velocity: 
delta_vO= [xdotO ; ydotO ; zdotO ] ;  %m/s 
%  Other  variables: 
t=t '  ; 

psi=n*t;  %  normalized  time 


%  I.  State  transition  matrices: 

zer=zeros (size (t) ) ; 
one=zer ; 
one ( : ) =1 ; 

Phi_rr= [4-3*cos (psi)  zer  zer ; 6* ( sin (psi ) -psi )  one  zer;zer  zer  cos(psi)] 
%Eq.  3.36 

Phi_rv= [(l/n)*sin(psi)  (2/n)*(l-cos(psi))  zer;(2/n)*(cos(psi)-l) 

(4/n) *sin (psi) - (3/n) *psi  zer;zer  zer  (1/n) *sin (psi) ] ;  %Eq.  3.37 

Phi_vr= [3*n*sin (psi)  zer  zer; 6*n* (cos (psi) -one)  zer  zer; zer  zer  - 
n*sin(psi)];  %Eq.  3.38 

Phi_vv= [cos (psi)  2*sin(psi)  zer; -2*sin (psi)  -3*one+4*cos (psi)  zer;zer 
zer  cos (psi)];  %Eq.  3.39 


%  II.  Solution : 

%Position  vector: 

delta_r_vector=Phi_rr*delta_rO+Phi_rv*delta_vO ;  %Eq.  3.34 
%Velocity  vector: 

delta_v_vector=Phi_vr*delta_rO+Phi_vv*delta_vO ;  %Eq.  3.35 

%Scalar  state  variables: 
x=delta_r_vector ( 1 : length (t) ) ; 
y=delta_r_vector (length (t) +1 : 2*length (t) ) ; 
z=delta_r_vector (2*length (t) +1 : 3*length (t) ) ; 
x_dot=delta_v_vector ( 1 : length (t) ) ; 
y_dot=delta_v_vector (length (t) +1 : 2*length (t) ) ; 
z_dot=delta_v_vector (2*length (t) +1 : 3*length (t) ) ; 
x=x '  ; 
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y=y 

z=z 


%  III.  ROE  printouts: 

a_e=2*sqrt  (  (3*x0+2*ydot0/n)  ^2+  (xdotO/n)  ''2) 
x_d=4  *x0+2*ydot0 / n 
y_d0=y0-2*xdot0/n- (3/2) *n*x_d*tO 
beta0=atan2 (xdotO, 3*n*x0+2*ydot0 ) 
z_max=sqrt (z0^2+ (zdotO/n) ^2) 
psi0=atan2 (zO, zdotO/n) 


Transforming  Virtual  Chief  Parameters  to  Cartesian  Initial  Conditions 

If  the  Virtual  Chief  Parameters  are  the  input  for  a  particular  case,  then  the  script 
VCM_parameter_tool.m  plots  the  Virtual  Chief  trajectory  and  prints  the  equivalent 
Cartesian  initial  conditions  as  outputs.  These  Cartesian  initial  conditions  can  then  be 
used  as  inputs  to  rel_sat_motion5.m. 

The  Algorithm. 

The  algorithm  is  as  follows: 

•  1.  Calculations.  Using  two-body  properties  from  Appendix  D,  calculate 
known  constants,  epoch  conditions,  and  arrays  of  chief  true  anomaly  and  time. 

•  n.  Parameterized  motion  model.  Calculate  the  solution,  c-frame 
coordinates  for  each  time  step,  using  the  results  of  Chapter  IV. 


•  ni.  Plots.  Generate  figures  similar  to  those  of  rel_sat_motion5.m  showing 
the  relative  coordinates  and  the  trajectory. 


•  rV.  Equivalent  Cartesian  initial  conditions.  Compute  and  print  Cartesian 
initial  conditions  equivalent  to  the  input  Virtual  Chief  parameters. 

The  Code. 

This  is  VCM_parameter_tool.m: 

%  Relative  orbit  plotter 

g, 

o 

%  This  is  a  stand-alone  script  that  plots  relative  satellite  motion 
using 
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%  the  six  new  virtual-chief  method  parameters  as  inputs.  The  script 
prints 

%  the  equivalent  Cartesian  initial  conditions,  which  can  then  be  used  as 
%  inputs  to  rel_sat_motion .m. 

clc;  clear  all;  close  all;  format  long  g 

g, 

0 


%  INPUT: 

g, 

0 


Al=  ;  %m  (gives  the  order  of  magnitude  of  the  y-motion) 

phil=  ;  %rad  (if  phi2=0,  gives  the  approximate  starting  angle  in  the 
%x-y  plane,  measured  clockwise  from  the  y-axis;  affects  "skewness"  and 
%"pointiness"  of  the  x-y  projection;  for  trajectories  with  out-of-plane 
%motion,  affects  the  orientation  of  the  3D  (near-planar ) traj ectory) 

A2=  ;  %m  (A2*e_C  gives  the  order  of  magnitude  of  the  x-motion;  A2 
%also  scales  the  x-  and  y-drift  when  phi2  is  not  zero  or  pi;  when  phi2 
%equals  zero  or  pi,  the  coordinates  of  the  "center"  of  the  near-ellipse 
are 

% ( 0 , A2 , 0 ) . ) 

phi2=  ;  %rad  (if  A2=0,  this  has  no  effect;  otherwise,  gives  an 
%indication  of  drift  rate;  dominant  drift  is  along  the  y-direction,  and 
is 

%negative  for  0<phi2<pi,  and  positive  for  pi<phi2<2 *pi ;  maximum  drift 
for 

%phi2=pi/2  and  3*pi/2;  that  includes  the  t*sin (t) -type  drift  in  the 
%x-di recti on) 

z_max=  ;  %m  (gives  the  amplitude  of  the  z-motion) 

psiO=  ;  %rad  (for  nonzero  z  max  and  phil=0,  gives  the  orientation  about 
the 

%z-axis  of  the  (near-planar)  3D  trajectory,  measured  clockwise  from  the 
%negative  x-axis) 

e_C=0.001  ;  %chief  orbit  eccentricity 
a_C=7778137  ;  %m  (chief  semi -major  axis) 

M_C0=0  ;  %rad  (chief  mean  anomaly  at  epoch) 

N=4  ;  % (Desired  number  of  orbits) 

ss=pi/360  ;  %rad  (Desired  step  size,  in  increments  of  chief  true 
anomaly) 

mu=3 . 98600441*10^14;  %  (m*3)  /  (s''2)  (gravitational  parameter) 
t0=0  ;  %s  (epoch  time) 


g, 

0 


%  I.  CALCULATIONS: 
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n_C=sqrt  (mu/ (a_C''3)  )  ;  %rad/s  (chief  mean  motion) 
t_p=tO-M_CO/n_C;  %s  (chief  time  of  perigee  passage) 

P_C=2*pi/n_C;  %s  (chief  period) 

flag=0;  %Newton-Raphson  iteration  to  find  epoch  conditions 

E0=M_C0;  %  first  guess 

El=10; 

while  flag==0, 

if  abs  (El-EO)  <10^'  (-9)  , 
f lag=l ; 
else  E0=E1; 

E1=E0- (E0-e_C*sin (EO) -M_C0) / (l-e_C*cos (EO) ) ; 

end 

end 

E_C0=E1;  %rad  (chief  eccentric  anomaly  at  epoch) 
if  E_C0<0 

E_C0=E_C0+2*pi; 

end 

if  M_C0==0 
E_C0=0; 

end 

s=sin  (E_C0  )  *sqrt  ( l-e_C''2  )  /  ( l-e_C*cos  (E_C0  )  )  ; 
c= (cos (E_C0) -e_C) / (l-e_C*cos (E_C0) ) ; 

nu_C0=atan2 (s, c) ;  %rad  (chief  true  anomaly  at  epoch) 
if  nu_C0<0 

nu_C0=nu_C0+2  *pi ; 

end 

nudot_C0=n_C*  (l-e_C'^2)  (-3/2)  *  (l+e_C*cos  (nu_C0)  )  ^2;  %rad/s  (chief 

angular  rate  at  epoch) 

nu_C=nu_C0 : ss :nu_C0+2*N*pi;  %rad  (chief  true  anomaly  array,  allowed  to 

take  on  values  greater  than  2*pi) 

s2=sin  (nu_C)  .  *sqrt  ( l-e_C''2  )  .  /  ( l+e_C*cos  (nu_C)  )  ; 

c2= (e_C+cos (nu_C) ) . / ( l+e_C*cos (nu_C) ) ; 

E_C=atan2 ( s2 , c2 ) ;  %rad  (chief  eccentric  anomaly  array) 
i=l;  %atan2  is  defined  from  -pi  to  pi;  I  need  zero  to  2*pi 
while  i<=length (E_C) 
if  E_C (i) <0 

E_C (i) =2*pi+E_C (i) ; 

end 
i=i+l ; 

end 

M_C=E_C-e_C*sin (E_C) ;  %rad  (chief  mean  anomaly  array) 
t=M_C . /n_C+t_p;  %s  (time  array) 

i2=l;  %I  don't  want  to  reset  time  to  zero  or  negative  after  each  orbit 
while  i2<=N 
i3=2; 

while  i3<=length (t) 
if  t(i3)<t(i3-l) 

t (i3) =t (i3) +P_C; 

end 

i3=i3+l ; 

end 

i2=i2+l; 

end 


g, 

0 


%  PARAMETERIZED  MOTION  MODEL: 


164 


x=Al*sin (nu_C+n_C*t_p-phil ) +A2*sin (nu_C- 

n_C*t+n_C*t_p+phi2 ) + (Al/3) *sin (nu_C-2*n_C*t+n_C*t_p+phil ) - 
(3/2) *A2*sin (phi2) *n_C* (t-tO) . *sin (nu_C-n_C*t+n_C*t_p) ; 
y=Al*cos (nu_C+n_C*t_p-phil ) +A2*cos (nu_C- 

n_C*t+n_C*t_p+phi2 ) + (Al/3) *cos (nu_C-2*n_C*t+n_C*t_p+phil ) - 
(3/2) *A2*sin (phi2) *n_C* (t-tO) . *cos (nu_C-n_C*t+n_C*t_p) ; 
z=z_max*sin (n_C* (t-tO) +psiO)  ; 


%  PLOTS: 

g, 

0 


subplot (3,1,1); 
xlabel ( ' t ' ) 
ylabel ( ' x ' ) 
title ( ' Relative 
subplot  (3,1,2); 
xlabel  (  ' t '  ) 
ylabel ( ' y ' ) 
ylabel  (  ' z  '  ) 
subplot (3,1,3); 
xlabel ( ' t ' ) 


plot (t, X, ' b ' ) 

coordinates  vs. 
plot (t, y, 'b' ) 

plot ( t , z , ' b ' ) 


time ' ) 


figure 

subplot (3,1,1);  plot (nu_C, x, ' b ' ) 
xlabel ( ' nu_C ' ) 
ylabel ( ' x ' ) 

title (' Relative  coordinates  vs.  chief  true  anomaly') 
subplot ( 3 , 1 , 2 ) ;  plot(nu_C,y) 
xlabel ( ' nu_C ' ) 
ylabel ( ' y ' ) 

subplot  ( 3 , 1 , 3 ) ;  plot(nu_C,z) 
xlabel ( ' nu_C ' ) 
ylabel  (  ' z  '  ) 


figure 

subplot (2 , 2 , 1 ) ;  plot (x, y, ' b ' ) 
axis  equal 
xlabel ( ' X ' ) 
ylabel ( ' y ' ) 

hold  on;  plot (x ( 1 ) , y ( 1 ) , ' ro ' ) 
plot (x (15) , y (15) , ' rx' ) 
title (' Relative  trajectory') 
subplot (2 , 2 , 2 ) ;  plot(x,z) 
axis  equal 
xlabel ( ' X ' ) 
ylabel  (  ' z  '  ) 

hold  on;  plot (x ( 1 ) , z  ( 1 ) ,  ' ro ' ) 
plot(x(15) ,z(15) , 'rx') 
subplot (2, 2, 3) ;  plot(y,z) 
axis  equal 
xlabel ( ' y ' ) 
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ylabel  (  ' z  '  ) 

hold  on;  plot (y ( 1 ) , z ( 1 ) , ' ro ' ) 

plot(y(15) ,z(15) , 'rx') 

subplot (2 , 2 , 4 ) ;  plot3(x,y,z) 

axis  equal 

grid  on 

xlabel ( ' X ' ) 

ylabel ( ' y ' ) 

zlabel  (  ' z  '  ) 

hold  on;  plot3 (x ( 1 ) , y ( 1 ) , z  ( 1 ) ,  '  ro  '  ) 
plot3 (x(15) ,y(15) , z (15) , 'rx' ) 


g, 

0 


%  EQUIVALENT  CARTESIAN  INITIAL  CONDITIONS: 

g, 

0 


xO=AI*sin (nu_CO+n_C*t_p-phil ) +A2*sin (nu  CO- 

n_C*tO+n  C*t  p+phi2) + (1/3) *AI*sin (nu  C0-2*n_C*t0+n_C*t  p+phil) 
yO=AI*cos (nu_CO+n_C*t_p-phil ) +A2*cos (nu  CO- 

n_C*tO+n^C*t  p+phi2 ) + ( I /3 ) *AI *cos (nu  C0-2*n_C*t0+n_C*t  p+phil) 
zO=z_max*sin (psiO) 

xdotO=-A2*n_C*sqrt  (1-  (3/4)  *  (sin  (phi 2)  )  ''2)  *cos  (nu  CO- 

n  C*tO+n_C*t  p+atan2 (-sin (phi2) , 2*cos (phi2) ) ) - (2/3) *Al*n_C*cos (nu_C0- 
2 *n_C*tO+n_C*t  p+phil) +nudot_CO*Al*cos (nu_C0+n_C*t  p- 
phil) +nudot_C0*A2*cos (nu_C0- 

n_C*tO+n_C*t  p+phi2) + (1/3) *nudot_CO*Al*cos (nu_C0-2*n_C*t0+n_C*t  p+phil) 
ydotO=A2*n_C*sqrt (1- (3/4) * (sin (phi 2) ) ^2) *sin (nu  CO- 

n_C*tO+n_C*t_p+atan2 (-sin (phi2) , 2*cos (phi2) ) ) + (2/3) *Al*n_C*sin (nu_C0- 
2 *n_C*tO+n_C*t  p+phil) -nudot_CO*Al*sin (nu_C0+n  C*t  p-phil)- 
nudot_C0*A2*sin (nu_C0-n  C*tO+n_C*t  p+phi2) - (1/3) *nudot_CO*Al*sin (nu_C0- 
2  *n_C*tO+n_C*t_p+phil ) 
zdotO=z_max*n_C*cos (psiO) 
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Appendix  G:  A  Decreasing  Function  of  Time 

This  section  will  prove  that  Am  (t)  -  is  always  a  decreasing  function  of  time 

for  values  of  less  than  about  0.3106016499352. 

From  Eq.  (4.44), 


Am  (t)  -  n^t  =  Vp  (t)  -  2n^t  + 

Thus,  our  goal  is  to  find  when  v  (t)  -  2n^t  +  is  a  decreasing  function  of  time. 

In  other  words,  to  find  when 


^(^(0-2v+Vp)<o 


(G.l) 


Evaluating, 


v(t)-2Mp  <  0 


Applying  the  property  of  two-body  orbits  contained  in  Eq.  (D.8)  and  rearranging, 

_3  /  2 

2^c>^c(l~^c  )  (0]) 


Assuming  an  elliptical  orbit  (i.e.,  less  than  unity). 
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>(l  +  ^’c  cos[vc(0]f 


(G.2) 


If  Inequality  (G.2)  is  true  when  the  right-hand  side  is  a  maximum,  then  it  is  always 
true.  The  right-hand  side  is  a  maximum  at  chief  perigee,  when  cos[v^  (t)]  equals  unity. 
Thus,  we  must  determine  for  what  values  of  the  following  is  true: 

2(l-e,‘f>(l  +  e,f  (G3) 

Of  course,  is  already  confined  to  the  domain  0<e^  <1.  As  increases,  the 
left-hand  side  of  Inequality  (G.3)  decreases,  and  the  right-hand  side  increases.  Thus,  the 
maximum  value  of  for  which  Inequality  (G.l)  is  always  true  can  be  found  from 

After  squaring,  expanding,  and  collecting  terms,  we  have 

+  -  -1-3  =  0 

The  above  polynomial  has  only  one  real  root,  0.310601649935225.  (The  complex 
roots  are  1.34469917503239  ±  0.778750642891565/ .) 

Thus,  for  less  than  or  equal  to  0.310601649935225,  Au[t)-n^t  is  a 
decreasing  function  of  time. 
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